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ABSTRACT 

We use cosmological simulations to study the effects of self-interacting dark matter (SIDM) 
c"| on the density profiles and substructure counts of dark matter halos from the scales of spi- 

t's [ ral galaxies to galaxy clusters, focusing explicitly on models with cross sections over dark 

matter particle mass cr/m = 1 and 0.1 cm 2 /g. Our simulations rely on a new SIDM N-body 
algorithm that is derived self-consistently from the Boltzmann equation and that reproduces 
analytic expectations in controlled numerical experiments. We find that well-resolved SIDM 
halos have constant-density cores, with significantly lower central densities than their CDM 
counterparts. In contrast, the subhalo content of SIDM halos is only modestly reduced com- 
pared to CDM, with the suppression greatest for large hosts and small halo-centric distances. 
Moreover, the large-scale clustering and halo circular velocity functions in SIDM are effec- 
tively identical to CDM, meaning that all of the large-scale successes of CDM are equally 
well matched by SIDM. From our largest cross section runs we are able to extract scaling 
relations for core sizes and central densities over a range of halo sizes and find a strong cor- 
relation between the core radius of an SIDM halo and the NFW scale radius of its CDM 
counterpart. We construct a simple analytic model, based on CDM scaling relations, that cap- 
tures all aspects of the scaling relations for SIDM halos. Our results show that halo core 
densities in cr/m = 1 cm 2 /g models are too low to match observations of galaxy clusters, 
low surface brightness spirals (LSBs), and dwarf spheroidal galaxies. However, SIDM with 
cr/m ~ 0.1 cm 2 /g appears capable of reproducing reported core sizes and central densities of 
dwarfs, LSBs, and galaxy clusters without the need for velocity dependence. Higher resolu- 
tion simulations over a wider range of masses will be required to confirm this expectation. We 
discuss constraints arising from the Bullet cluster observations, measurements of dark matter 
density on small-scales and subhalo survival requirements, and show that SIDM models with 
cr/m ~ 0.1 cm 2 /g ~ 0.2barn/GeV are consistent with all observational constraints. 
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1 INTRODUCTION 

There is significant evidence that some form of dark matter dom- 
inates the gravitating ma ss in the universe an d its abundance is 
known to great precision dKomatsu et alj|2oTlh . The most popular 
candidate for dark matter is the class of weakly interacting massive 
parti cles (WIMPs), of which supersymmetric neutralinos are exam - 
ples dSteigman & Turnenll985l ; lGriestlll988l;|jungman et alj|l996h . 
WIMPs are stable, with negligible self-interactions, and are non- 
relativistic at decoupling ("cold"). It is important to recognize that 
of these characteristics, it is primarily their coldness that is well 
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tested via its association with significant small-scale power. Indeed, 
WIMPs are the canonical Cold Dark Matter (CDM) candidate. Cos- 
mological models based on CDM repr oduce the spatia l clustering 
of galaxies on large scales quite well teeid et al.ll2oToh and even 



the clustering of galaxies on < 

expected for CDM subhalos jKravtsov et alj|2004l: IConrov et al.l 



1 Mpc scales appear s to match that 



l2006l : lTruiillo-Gomez et alj|201 ll:lReddick et alj201 



Beyond the fact that the universe appears to behave as ex- 
pected for CDM on large scales, we have few constraints on the 
microphysical parameters of the dark matter, especially those that 
would manifest themselves at the high densities associated with 
cores of galaxy halos. It is worth asking what (if anything) about 
vanilla CDM can change without violating observational bounds. 
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In this paper we use cosmological simulations to explore the ob- 
servational consequences of a CDM particle that is strongly self- 
interacting, focusing specifically on the limiting case of velocity- 
independent, elastic scattering. 

Dark matter particles with appreciable self-interactions have 
been discussed in t he literature for mor e than two decades 
dCarlson et al] 1 19921. iMachacek et all 1 19931: Ide Laix et al.l 1 19951 : 
Spergel & Steinh ardt 2000; Firm ani et al Jl 20001) . and are now rec- 
ognized as generi c consequences of hidden-sector extensions to the 
Standard Model dPospelov et al]|2008l; lArkani-Hamed et alj|2009l; 
Ackerman etail 120091 ; iFeng et alj l2009t l20ld : lLoeb & Weineil 



201 lh . Importantly, even if dark sector particles have no couplings 



to Standard Model particles they might experience strong in terac- 
tions wit h themselves , mediated by dark gauge bosons (see iFend 
l20ld and |Peteill2012l for reviews). The implication is that astro- 
physical constraints associated with the small-scale clustering of 
dark matter may be the only way to test these scenarios. 

Phenomenologically, self-interacting dark matter (SIDM) is 
attractive because it offers a means to lower the central den- 
sities of galaxies without destroying the successes of CDM 
on large scales. Cosmological simulations that contain only 
CDM indicate that dark-matter halos should be cuspy and with 
(high) c oncentrations that correlate with the collapse time of 
the h alo dNavarro et aflll99l iBuTlock et al-lEoOll : IWechsler et all 
2002). This is inconsistent with observations of galaxy rotation 
curves, which show that galaxies are less conc entrated and less 



cuspy than predicted in C D M simulations (e.g. iFlores & 



1994 : ISimon et al] 120051: iKuzio de Narav etail 120081: 



20ld: button et al] l201ll: IKuzio de Narav & Spekkem 



de Blokl 



Oh et al] 1201 ll: IWalker & Penarrubial 1201 it ISalucci et 



k 



2011 



201/ 



Castignani et "atl l2012h . Even for clusters of galaxies, the density 



profiles of the host dark-matter halos appear in a number of cases 
to be shallower than predicted by CDM-only structure simula- 
tions, with the total (dark matter + baryons) density profile in a 
closer match to the CDM prediction for the dark matter alone (e.g . 
Sand et alj20oll2008l : lNewman et alj2009ll201 lhlCoe et alj2012t 
Umetsu et al]|2012h . 

One possible answer is feedback. In principle, the expul- 
sion of gas from galaxies can result in lower dark matter den- 
sities compared to dissipationless simul ations, and thus bring 



CDM m odels i n line with observatio ns (Governato et al 



Oh et al] I20TI 



Governato et al] 



Pontzen & Governatd l201ll : iBrooketal 



2010; 



2012 



2012h . However, a new level of concern ex- 



ists for dwarf spheroidal galaxies dBovlan-Kolchin et al] l2011al ; 
iFerrero et aD 1201 it iBovlan-Kolchin et al] 1201 lbl) . Systems with 
M* ~ 10 6 M Q appear to be missing ~5x 10 7 M fl of dark mat- 
ter com pared to standard CDM expectations (Bovlan-Kolchi net al] 
1201 lbh . It is difficult to understand how feedback from such a tiny 
amount of star formation could have possibly blown out enough 
gas to reduce the densities of dwarf spheroidal galaxies to the 
level required to match observations dBovlan-Kolchin et al.ll2011b ; 



2012 



iTevssier et all 120121 : IZolotov et al] |2012| ; IPenarrubia et ~ 
Garris on- Kimmel et al., in preparation). 

Spergel & Steinhar dt(2000) were the first to discuss SIDM in 
the co ntext of the central density problem (see also iFirmani et al.l 
2000). The centers of SIDM halos are expected to have constant 
density isothermal cores that aris e as kinetic energy is transmit- 
ted from the hot outer halo inward (Balbera et al. 2002: IColm et al] 
l2002tlAhn & Shapirdl2005l ; lKoda & Shapirdl201 lh . This can hap- 
pen if the cross section over mass of the dark matter particle, 
a/m, is large enough for there to be a relatively high probabil- 
ity of scattering over a time t agc comparable to the age of the halo: 



vary with local dark matter density p(r) as a function of radius r in 
a dark halo as 



T(r) ~ p(r)(a/m)v rlaa (r) , 



(1) 



up to order unity factors, where v rms is the rm s speed of dark-matter 
particl es. Based on rough analytic arguments, Spergel & Steinhardtl 
(2000) suggested a /m ~ 0.1 — 100 cm 2 /g would produce observ- 
able consequences in the cores of halos. 

Numerical simulations have c onfirmed the expected 
phenomenology o f core formation dBurker I I2000I) though 
iKochanek & White! d200d) emphasized the possibility that SIDM 
halos could eventually become more dense than their CDM 
counterparts as a result of eventual heat flux from the inside 
out (much like core collapse globular clusters). However this 
process is suppressed when mer ging from hierarchic al formation 
is included (for a discussion see I Ann & Sha piro 2005). We do not 
see clear signatures of core collapse in the halos we analyzed for 
a/m = 1 cm 2 /g. 

The first cosmological simulat ions aimed at un derstanding 
dwarf densities were performed by iDave et al] d200lh who used 
a small volume (AhT 1 Mpc on a side) in order to focus com- 
putational power on dwarf halos. They concluded that a/m = 
0.1 — 10cm 2 /g came close to reproducing core densities of small 
galaxies, favoring the upper end of that range but not being able 
to rule out the lowe r end due to resolution. Almost concurrently, 
lYoshida etall {2000 ) ran cosmological simulations focusing on the 
cluster-mass regime. Based on the estimated core size of cluster 
CL 0024+1654, they concluded that cross sections no larger than 
~ 0.1 cm 2 /g were allowed, raising doubts that constant-cross- 
section SIDM models could be consistent with observations of both 
dwarf galaxies and clusters. 

These concerns were echoed by iMiralda-Esc ude ( 2002) who 
suggested that SIDM halos would be signific antly more spheri- 
cal tha n observed for galaxy clusters. Similarlv. lGnedin & Ostrikeil 
feooih argued that SIDM would lead to excessive sub halo evapo- 
ration in galaxy clusters. More recently, the merging cluster system 
known as the Bullet Cluster has been used to de rive the limits (68% 
CL.) a/m < 0.7cm 2 /g dRandall et al1l2008l) based on evapora- 
tion of dark matter from the subcluster and o/m < 1.25cm 2 /g 
dRandall et al]|2008l) based on the observed lack of offset between 
the bullet subcluster mass peak and the galaxy light centroid. In 
order to relax this apparent tension between what was required 
to match dwarf densities and the observed properties of galaxy 
clusters, velocity dependent cross sections that diminish the ef- 
fects of self-interaction in cluster environments have been con- 



sidered (Firmani et al. 200d ; IColm et al ]|2002l ; lFeng etal]|2009l: 
lLoeb & Weineill201 ltlVogelsberger et alj2012h . 

There are a few new developments that motivate us to revisit 
constant SIDM cross sections on the order of a/m ~ 1 cm 2 /g. For 
example, the cluster (CL 0024+1654) used bv lYoshida et al] 1 2000) 
to place one of the tightest limits at a/m — 0. 1, is now recognized 
as an ongoing merger along the line of s ight ( Czoske et al .11200 it 
20021 ; IZhang et aljfeoOSt Ijee et al.l [20071 : |jeell20ld : lUmetsu et alj 
2010). This calls into question its usefulness as a comparison case 
for non-merging cluster simulations. In a companion paper (Peter, 
Rocha, Bullock and Kaplinghat, 2012) we use the same simulations 
described here to show that published constraints on SIDM based 
on halo shape comparisons are significantly weaker than previ- 
ously believed. Further, the results presented below clearly demon- 
strate that the tendency fo r subhalos to evaporate in SIDM models 
dGnedin & Ostrikeill200ll) is not significant for cr/m ~ 1 cm 2 /g. 
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Finally (and related to the previous point), the best numerical anal- 
ysis of the Bullet Cluster dRandall et alj|2008l) used initial cluster 
density profiles that were unmotivated cosmologically with central 
densities about a factor of two too high for the SIDM cross sec- 
tions considered (producing a scattering rate that is inconsistently 
high). Based on this observation, the bullet cluster constraint based 
on evaporation of dark matter from the subcluster should be re- 
laxed since the amount of subcluster mass that becomes unbound 
is directly proportional to the density of dark matter encountered in 
its orbit. Moreover, their model galaxies were placed in the clus- 
ter halo potentials without subhalos surrounding them, an assump- 
tion (based on analytic estimates for SIDM subhalo evaporation) 
that is not supported by our simulations. This could affect the con- 
straints based on the (lack of) offset between dynamical mass and 
light. Thus we believe that the bullet cluster constraints as discussed 
above are likely only relevant for models with a/m > lcm 2 /g. 
However, the constraints could be made significantly stronger by 
comparing SIDM predictions to the densities inferred from the con- 
vergence maps since the central halo densities for a/m ~ 1 cm 2 /g 
are significantly lower than the CDM predictions, as we show later. 

Given these motivations, we perform a set of cosmologi- 
cal simulations with both CDM and SIDM. For SIDM we ran 
<j/m = l and 0.1 cm 2 /g models (hereafter SIDMi and SIDMo.i), 
i.e., models that we have argued pass the Bullet cluster tests. Our 
simulations provide us with a sample of halos tha t span a mass 
range much larger than either iDave et al or lYoshida et al.l 

(2000) both with and without self-interactions. 

One of the key findings from our simulations is that the core 
sizes are expected to scale approximately as a fixed fraction of the 
NFW scale radius the halo would have in the absence of scatter- 
ings. We can see where this scaling arises from a quick look at 
Equation [T] This equation allows us to argue that the radius (ri) 
below which we expect dark matter particles (on average) to have 
scattered once or more is set by: 

vl x 

Ps/(r/r a )w rms oc -g— /(ri/r s )Kiax = constant , (2) 

where f(x) is the functional form of the NFW (or a related) den- 
sity profile. In writing the above equation we have assumed that the 
density profile for SIDM is not significantly different from CDM at 
ri, something that we verify through our simulations. Now, since 
CDM enforces a V ma x — r max relation such that V ma x oc r]^~ , 
we see that the solution to ri/r s is going to be only mildly depen- 
dent on the halo properties. We develop an analytic model based 
on this insight later, but this is the underlying reason for why we 
find core sizes to be a fixed fraction of the NFW scale radius of the 
same halo in the absence of scatterings. 

The major conclusion we reach based on the simulations and 
the analytic model presented here is that a self-interacting dark 
matter model with a cross-section over dark matter particle mass 
~ 0.1 cm 2 /g would be capable of reproducing the core sizes 
and central densities observed in dark matter halos at all scales, 
from clusters to dwarf spheroidals, without the need for velocity- 
dependence in the cross-section. 

In the next section, we discuss our new algorithm to compute 
the self-interaction probability for N-body particles, derived self- 
consistently from the Boltzmann equation. We discuss this new 
algorithm in detail in Appendix [A] In £|2] we show how this al- 
gorithm is impl emented in the publicly available code GADGET-2 
(Springel 2005). We run tests that show that our algorithm gets the 
correct interaction rate and post-scattering kinematics. The results 
of these tests are in ij3] The cosmological simulations with this new 



algorithm are described in detail in jj4] In £15.11 we provide some 
preliminary illustrations of our simulation snapshots and in !5.2l we 
demonstrate that the large-scale statistical properties of SIDM are 
identical to CDM. In ^5 . 3 1 we present the properties of individual 
SIDMi and SIDMo.i halos and compare them to the their CDM 
counterparts. In i]5.4| we discuss the subhalo mass functions in our 
SIDM and CDM simulations and show that SIDMi subhalo mass 
functions are very close to that of CDM in the range of halo masses 
we can resolve. We provide scaling relations for the SIDMi halo 
properties in ij6]and in i|7]we present an analytic model that repro- 
duces these scaling relations as well as the absolute densities and 
core radii of SIDMi halos. We use these scaling relations and the 
analytic model to make a broad-brush comparison to observed data 
in Sj8] We present a summary together with our final conclusions in 

m 



2 SIMULATING DARK MATTER SELF INTERACTIONS 

Our simulations rely on a new algorithm for modeling self- 
interacting dark matter with N-body simulations. Here we intro- 
duce our approach and provide a brief summary. In Appendix lAlwe 
derive the algorithm explicitly starting with the Botlzmann equa- 
tion and give details for general implementation. 

In N-body simulations, the simulated (macro)particles repre- 
sent an ensemble of many dark-matter particles. Each simulation 
particle of mass m p can be thought of as a patch of dark-matter 
phase-space density. In our treatment of dark matter self-scattering, 
the phase space patch of each particle is represented by a delta func- 
tion in velocity and a spatially extended kernel W(r, h s i), smooth- 
ing out the phase space in configuration space on a self-interaction 
smoothing length h sl - The value of h s i needs to be set by consid- 
ering the physical conditions of the problem (see ^3) as it specifies 
the range over which N-body particles can affect each other via 
self-interactions. In principle, h s i could be different for each parti- 
cle and vary depending on the local density, but in the simulations 
presented here we fix h s i to be the same for all particles in a given 
simulation, setting the size of hsi according the lowest densities at 
which self-interactions are effective for a given cross section. 

When two phase-space patches overlap, we need to calculate 
the pairwise interaction rate between them. We do so by consid- 
ering the "scattering out" part of the Boltzmann collision term in 
Equation iAH and Eqs. l |A81 >-( lAT3~b . The implied rate of scattering 
of an N-body particle j off of a target particle i of mass m p is 

r (*|j) = (cr/m)m p \vi - Vj\gji , (3) 

where gji is a number density factor that accounts for the overlap 
of the two particles' smoothing kernels: 

g it = [ " d 3 x'W(\x'\,h Bi )W(\S Xjl + x'|,M . (4) 
Jo 

The probability that such an interaction occurs in a time step 8t is 

P(i\j)=T(i\j)5t, (5) 

and the total probability of interaction between N-body particles i 
and j is 

Ph = Em+mii. (6) 



1 This equation applies only if /i s ; is the same for both particles. See Ap- 
pendix A for the general form. 
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Specifically, Pij is the probability for a macroparticle representing 
a patch of phase space around (xj,Vj) to interact with a target 
particle representing a patch of phase space around (x;,Vi) in a 
time St. 

We determine if particles interact by drawing a random num- 
ber for each pair of particles that are close enough for the proba- 
bility of interaction to be greater than zero. If a pair does scatter, 
we do a Monte Carlo for the new velocity directions, populating 
these parts of the phase-space and deleting the two particles at their 
initial phase-space locations. Note that by virtue of populating the 
new phase space regions, we are taking care of the "scattering in" 
term of the collision integral in Equation JAlb . We avoid double 
counting by only accounting for Pij — Pji once during a given 
time-step St. In the limit of a large number of macroparticles, the 
total interaction probability for each particle i should approach 



Pi — y t Pij 



(7) 



We show in §3 that this approach correctly reproduces the expected 
number of scatterings in a idealized test case. 

Our method for simulating scattering differs from previous ap- 
roach es in a few key ways. It is most similar to that of lDave et al.l 
200 ll) in that we both directly consider interactions between pairs 
of phase-space patches and rely on a scattering rate similar in form 
to Equation[3] The difference is that their geometric factor gji is not 
the same — our factor arises explicitly from the overlap in patches of 
phase space between neighboring macroparticles, as derived from 
the collision term in the Boltzmann equation (see Appendix A for 
details). Other authors determine the scattering rate T of individ- 
ual phase-space patches based on estimates of the local mass den- 
sity (typically using some number of nearest neighbors or using 
an SPH kernel). The Monte Carlo is then based on an estimated 
scattering rate of an individual particle on the background, and a 
scattering partner is onl y chosen after a scattering event is deter- 
mined to have occurred (Kochan ek et al.l200 0; Yoshida ~et alj2000l ; 
ICoh'n et alj2002MRandall et alj200l) . The scattering probability in 
this latter approach is not symmetric. For macroparticles of iden- 
tical mass, P(i\j) = P(j\i) explicitly in our approach, but not 
the other approach because the density estimated at the position of 
macroparticle i need not be the same as that estimated at the posi- 
tion of particle j. In the future, there should be a direct comparison 
among these scattering algorithms to determine if they yield con- 
sistent results. 

We have implemented our algorithm in the publicly available 
versio n of the cosmological simulation code GADGET-2 dSpringell 
l2005t) . GADGET-2 computes the short-range gravitational interac- 
tions by means of a hierarchical multipole expansion, also known 
as a tree algorithm. Particles are grouped hierarchically by a re- 
peated subdivision of space, so their gravitational contribution can 
be accounted by means of a single multipole force computation. A 
cubical root node encompasses the full mass distribution. The node 
is repeatedly subdivided into eight daughter nodes of half the side 
length each (an oct-tree) until one ends up with "leaf" nodes con- 
taining single particles. Forces for a given particle are then obtained 
by "walking" the tree, opening nodes that are too close for their 
multipole expansion to be a correct approximation to their gravita- 
tional contribution. In GADGET-2, spurious strong close encoun- 
ters by particles are avoided by convolving the single point particle 
density distribution with a normalized spline kernel ("gravitational 
softening"). 

To implement our algorithm, we take advantage of the tree- 
walk already build in GADGET-2, computing self interactions dur- 
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Figure 1. Fraction of the expected total number of interactions that are com- 
puted in our test simulation as a function of the self-interaction smoothing 
length. The self-interaction cross section for each run is shown in units of 
cm 2 /g in the legend. The code converges to the expected number of interac- 
tions when the smoothing length approaches the background inter-particle 
separation, i.e. when /i s i(pbg/"^p) 1,/3 ii 0.2. 



ing the calculation of the gravitational interactions. For this to work 
we have to modify the opening criterion such that nodes are opened 
if they are able to have particles closer than 2ft s i from a target scat- 
terer (or hi + h j if particles have different self-interaction smooth- 
ing lengths). When computing the probabilit y of interaction we use 
the sa me spline kernel used in GADGET-2 ( Monaghan & Lattanziol 
Il985h . defined as 



1 - > 



+ 6| 



W { r,h) = — 3 



2(1-*) > 
0, 



< - < - 

u — h — 2 ' 

< 1, (8) 



h 

> 1. 



If a pair interacts we give both particles a kick consistent with 
an elastic scattering that is isotropic in the center of mass frame. 
The post-scatter particle velocities are 



v = v c + 



mi 



mo + mi 
mo 

mo + mi 



-Ve, 



Ve, 



(9) 



where v c is the center of mass velocity, V is the relative speed of 
the particles (conserved for elastic collisions) and e is a randomly 
chosen direction. 

The time-step criterion is also modified to assure that the scat- 
tering probability for any pair of particles is small, P — F St « 
1. An individual particle time-step is decreased by a factor of 2 if 
during the last tree-walk the maximum probability of interaction for 
any pair involving such a particle was P ma x > 0.2. Once a particle 
time-step is modified due to the previous restriction, if P max < 0.1 
for such a particle and its current time-step is smaller than the one 
given by the standard criterion on GADGET-2, we increase it by a 
factor of 2. 
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Table 1: Simulations discussed in this paper. 
Name Volume Number of Particles Particle Mass Force Softening Smoothing Length Cross-section 

£box [ft -1 Mpc] N p m p [/i -1 M0] e[/i _1 kpc] h si [h' 1 kpc] cr/m [cm 2 /g] 



CDM-50 


50 


512 3 


6.88 x 10 7 


1.0 







CDM-25 


25 


512 3 


8.59 x 10 6 


0.4 







CDM-Z11 


(3i?vir)* 


2.5 x 10 6 * 


1.07 x 10 6 * 


0.3 







CDM-Z12 


(3i? vir )* 


5.6 x 10 7 * 


1.34 x 10 5 * 


0.1 







SIDMo.i-50 


50 


512 3 


6.88 x 10 7 


1.0 


2.8 e 


0.1 


SIDMo.i-25 


25 


512 3 


8.59 x 10 6 


0.4 


2.8 e 


0.1 


SIDM .i-Zll 


(3-Rvlt)* 


2.5 x 10 6 * 


1.07 x 10 6 * 


0.3 


2.8 e 


0.1 


SIDMo.i-Z12 


(3R vi r)* 


5.6 x 10 7 * 


1.34 x 10 5 * 


0.1 


1.4 e 


0.1 


SIDMi-50 


50 


512 3 


6.88 x 10 7 


1.0 


2.8 e 


1 


SIDMi-25 


25 


512 3 


8.59 x 10 6 


0.1 


2.8 e 


1 


SIDMi-Zll 


(3-Rvir)* 


2.5 x 10 6 * 


1.07 x 10 6 * 


0.3 


2.8 e 


1 


SIDMi-Z12 


(3R vi r)* 


5.6 x 10 7 * 


1.34 x 10 5 * 


0.1 


1.4 e 


1 



*Note: The Zl 1 and Z12 runs are zoom simulations with multiple particle species concentrating on halos of mass M v ; r = 5 x 10 11 Mq and 1.0 X 10 12 
Mq, respectively (no h). The volumes listed refer to the number of virial radii used to find the Lagrangian volumes associated with the zoom. The particle 

properties listed are for the highest resolution particles only. 



3 TEST OF THE SIDM IMPLEMENTATION 

Before performing cosmological simulations, we carried out a con- 
trolled test of the implementation in order to make sure the scatter- 
ing rate and kinematics are correctly followed in the code, and to 
determine the optimum value of the SIDM softening kernel length 
h B i. The simplest and cleanest scenario for testing our implemen- 
tation consists of a uniform sphere of particles moving through a 
uniform field of stationary background particles. The coordinate 
system is defined such as the sphere is moving along the positive 
z-direction with constant velocity v s . The particles forming the 
sphere and the particles forming the background field are tagged 
as different types within the code and here we will refer to them 
simply as sphere (s) and background (bg) particles respectively. 
We only allow scatterings involving two different types of particles 
(i.e. sphere-background interactions only) and turn off gravitational 
forces among all of the particles. Furthermore all particles have the 
same mass m p . 

The expected number of interactions for this case is given by 



N exp (t) = Pij = N s {a/m)p hg v s t (10) 

where N s is the total number of Sphere particles, pb g is the density 
of the background field and t is the elapsed time from the begin- 
ing of the simulation. From this experiment we have found that the 
number of interactions computed by the code depends on the self- 
interaction smoothing length h B1 (see Figure [T), which is fixed to 
be the same for all particles in this test. The number of interac- 
tions converges to the expected value given by Equation < | 1 0b as /i s i 
becomes comparable to the background inter-particle separation, 
specifically when /isiCpbg/fip) 1 ' 73 > 0.2. For h s i(pi, s /m p ) 1 ^ 3 > 
0.5 the accuracy of the algorithm does not improve by much and 
the time of the calculations increases rapidly, oc ft 3 ;. Apart from 
the expense, using larger values of h s i would lead to increasingly 
non-local interactions among particles, which is inconsistent with 
the model under consideration. 

We also check the kinematics of the scatters in this test simula- 
tion and describe the results in Appendix|B] The resulting kinemat- 
ics and number of interactions from our test simulation agrees well 



with the expectations from the theory as long as h a i (pb g /m p ) 1//3 > 
0.2. 



4 OVERVIEW OF COSMOLOGICAL SIMULATIONS 

We initialize our cosmological sim ulations using the M ulti-Scale 
Initial Conditions (MUSIC) code of lHahn & Abel l bOllI) . We have 
a total of four initial condition sets, each run with both CDM 
and SIDM. The first two are cubic volumes of 25/i _1 Mpc and 
50/i -1 Mpc on a side, each with 512 3 particles. As discussed be- 
low, these simulations allow us to resolve the structure of a statisti- 
cal sample of group (~ 10 13 M ) and cluster (~ 10 14 Mq) halos. 

The second two initia l conditions concentr ate computational 
power on zoom regions dKatz& White] 1 19931) drawn from the 
50h~ Mpc box, specifically aimed at exploring the density struc- 
ture of two smaller halos, one with virial massQ Af v ir = 7.1 x 
lO 11 /^ 1 M Q = 1 x 10 12 M (Z12) and one with M vir = 3.5 x 
1O u /i _1 M = 5 x 10 11 Mq (Zll). The Z12 run in particular is 
fairly high resolution, with more than five million particles in the 
virial radius. Table 1 summarizes the simulation parameters. The 
cosmology used is based on WMAP7 results for a ACDM Uni- 
verse: h = 0.71, fi m = 0.266, Oa = 0.7 34, D. h = 0.0449, 
n s = 0.963, a 8 = 0.801 dKomatsu et alj201ll) . 

Each of our four initial conditions has been evolved from red- 
shift z = 250 to redshift z — with collisionless dark matter 
(labeled CDM) and with two types of self-interacting dark mat- 
ter: one with cr/m = 1 cm 2 /g (labeled SIDMi) and another 
with cr/m = 0.1 cm 2 /g (labeled SIDM0.1). We can use the 
same initial conditions for CDM and SIDM because at high red- 
shift the low densities and low relative velocities of the dark mat- 
ter make self-interactions insignificant. Table [2] list all the simu- 
lations used for this study and detail their force, mass, and self- 
interaction resolution. In addition to the simulations listed in the 
table, we also ran the cosmological boxes with SIDM cross sec- 
tions cr/m — 0.03 cm 2 /g. We do not present results from these low 



2 We define M v i r as M v i r = |7rp(,A v i r (.z)r 3 ir , and r v i r as p(r v i r ) = 
A v i r (,z)p(,. Where p(r v - lr ) denotes the overdensity within r v i r , is the 
background density and A v ; r the virial overdensity. 
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Figure 2. Top: Large scale structure in CDM (left) and SIDMi (right) shown as a 50x50 h _1 Mpc slice with 10 h _1 Mpc thickness through our cosmological 
simulations. Particles are colored according to their local phase-space density. There are no visible differences between the two cases. Bottom: Small scale 
structure in a Milky Way mass halo (Z12) simulated with CDM (left) and SIDMi (right), including all particles within 200/i _1 kpc of the halo centers. The 
magnitude of the central phase-space density is lower in SIDM because the physical density is lower and the velocity dispersion is higher. The core of the 
SIDM halo is also slightly rounder. Note that substructure content is quite similar except in the central regions 



cross section runs here because no core density differences were re- 
solved within the numerical convergence radii of our simulations. 

As shown in ij3]the self-interaction smoothing length h s i must 
be larger than 20% the inter-particle separation in order to achieve 
convergence on the interaction rate. All the work for this paper was 
done with a fixed h s i for all particles carefully chosen for each sim- 
ulation so that the self-interactions are well resolved at densities a 
few times to an order of magnitude lower than the lowest densities 
for which self-interactions are significant. We have run the cosmo- 
logical boxes with different choices for h s i (changes by factors of 
2 to 4) and have found that our results are unaffected. We have 
also run tests on isolated halos with varying smoothing lengths and 



again find that the effects of self-interactions are robust to reason- 
able changes in h s i . 



All of our halo catalogs and density profiles are derived 
using the publicly availab le code Amiga Halo Finder (AHF) 
jKnollmann & Kn ebe 2009). 
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Figure 3. Large-scale characteristics Left: Dark matter two-point correlation functions from our CDM-50 (CDM-25) and SIDMi-50 (SIDMi-25) simulations 
in black (grey) and blue (cyan) colors respectively. There are no noticeable difference between the CDM and SIDMi dark matter clustering over the scales 
plotted. Right: Cumulative number density of dark matter halos as a function of their maximum circular velocity (V max ) at different redshifts for our CDM-50 
(solid) and SIDMi-50 (dashed) simulations. There are no significant differences in the V max functions of CDM and SIDMi at any redshift. 



5 SIMULATION RESULTS 
5.1 Preliminary Illustrations 

Before presenting any quantitative comparisons between our CDM 
and SIDM runs, we provide some simulation renderings in order to 
help communicate the qualitative differences. 

The upper panels of Figure[2]show a large-scale comparison: 
two (50 x 50 x 10) h' 1 Mpc slices from the CDM-50 and SIDMj- 
50 boxes side-by-side at z — 0. The structures are color-coded by 
local phase-space density (oc p/«r ms ). It is evident that there are 
no observable differences in the large-scale characteristics of CDM 
and SIDMi. We discuss this result in more quantitative terms in 
S]5.2l but of course this is expected. The SIDM models we explore 
do not have appreciable rates of interaction for densities outside the 
cores of dark matter halos. The upper panels of Figure f2] provide a 
visual reminder that the SIDM models we consider are effectively 
identical to CDM on larges scales. 

The differences between CDM and SIDM become apparent 
only when one considers the internal structure of individual ha- 
los. The lower panels of Figure [2] provide side-by-side images of a 
Milky- Way mass halo (Z12) simulated with CDM (left) and SIDMi 
(right). SIDM tends to make the cores of halos less dense and ki- 
netically hotter (see $5.3\ and these two differences are enhanced 
multiplicatively in the phase-space density renderings. The central 
regions of the host halo are also slightly rounder in the SIDM case 
(Peter et al. 2012). Importantly, the difference in substructure char- 
acteristics are minimal, especially at larger radii. We return to a 
quantitative description of substructure differences in i]5.4| 



5.2 Large Scale Structure and Halo Abundances 

Figure [3] provides a quantitative comparison of both the clustering 
properties (left) and halo abundance evolution (right) between our 
full-box CDM and SIDMi simulations. The left panel shows the 
two-point function of dark matter particles in both cosmological 
runs for CDM and SIDMi. There are no discernible differences 



between SIDM and CDM over the scales plotted, though of course 
the different box sizes (and associated resolutions) mean that the 
boxes themselves only overlap for a limited range of scales. For 
a given set of initial conditions, however, SIDM and CDM give 
identical results. 

The right panel of Figure [3] shows the cumulative number 
density of dark-matter halos (including subhalos) as a function of 
their peak circular velocity (Vmax) for the CDM-50 (solid) and 
SIDMi-50 (dashed) simulations at various redshifts. Remarkably, 
this comparison shows no significant difference either - indicat- 
ing that SIDM with cross sections as large as lcm 2 /g does not 
strongly affect the maximum circular velocities of individual halos. 
The two panels of Figure [3] demonstrate that for large-scale com- 
parisons, including analyses involving field halo mass functions, 
SIDM and CDM yield identical results. The implication is that ob- 
servations of large-scale structure are just as much a "verification" 
of SIDM as they are of CDM. 

5.3 Halo Structure 

Before presenting statistics on halo structure, we focus on six well 
resolved halos that span our full mass range M V1I — 5 x 10 11 — 
2 x 10 14 Mq, selected from our full simulation suite, including 
our two zoom-simulation halos (Z12 and Zll). Figures [4] through 
[6] show radial profiles for the density, circular velocity and velocity 
dispersion for all three dark matter cases. In each figure, black cir- 
cles correspond to CDM, green triangles to SIDMo.i, and blue stars 
to SIDMi. All profiles are shown down to the innermost resolved 
radius for which the average t wo-body relaxatio n time roughly 
matches the age of the Universe jPower et alj2003l) . 

We begin with the density profiles of halos shown in the six- 
panel Figure|4] For each ha lo in the CDM run we have fit an NFW 
profile ^Navarro et al J 19971) to its radial density structure: 

p 

PnfwO) = . " * . , (11) 
r(r s + r-y 

and recorded its corresponding scale radius r s . The CDM-fit r 3 
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value for each halo is given in its associated panel along with the 
halo virial mass. The radial profiles for each halo (in both the CDM 
and SIDM runs) are normalized with respect to the CDM r s value 
in the plot. This allows our full range of halo masses to be plotted 
on identical axes. 

The SIDM versions of each halo show remarkable similar- 
ity to their CDM counterparts at large radii. However, the SIDMi 
cases clearly begin to roll towards constant-density cores at small 
radii. The best resolved halos in the SIDMo.i runs also demonstrate 
lower central densities compared to CDM, though the differences 
are at the factor of ~ 2 level even in our best resolved systems. 
Clearly, higher resolution simulations will be required in order to 
fully quantify the expected differences between CDM and SIDM 
for a/m ~ 0.1 cm 2 /g. 

Fo r the SIDMi ca ses we can quantify the halo cores by fitting 
them to lBurkertl t l995) profiles 



pB(r) 



Phr'h 



(r + r h )(r 2 + rg) ' 



(12) 



These Burkert fits are shown as blue dashed lines. They are good 
fits for radii within r ~ 2 — 3 r s , but the quality of the fits gets 
worse at large radii. The blue arrows in each panel show the value 
of the best-fit Burkert core radius for the SIDMi halos. Note that 
the values are remarkably stable in proportion to the CDM r s value 
at r b ~ 0.7 r 3 . 

As explained in the fact that the SIDMi profiles are rea- 
sonably well characterized by a single scale-radius Burkert profile 
may be a lucky accident, only valid for cross sections near 1 cm 2 /g. 
It just so happens that for this cross section the radius where dark 
matter particles experience significant scattering sets in at r ~ r s 
(see Figure [7] and related discussion). For a smaller cross section 
(with a correspondingly smaller core) a multiple parameter fit may 
be necessary. Given the beginnings of very small cores we are see- 
ing in the SIDMo.i runs, it would appear that we would need one 
scale radius to define an r s bend and a second scale radii to define 
a distinct core. 

Another qualitative fact worth noting is that the density pro- 
files of the SIDMi halos overshoot the CDM density profiles near 
the Burkert core radius (not as much as the Burkert fits do, but 
the difference in the data points is noticeable). This is due to the 
fact that as particles scatter in the center, those that gain energy 
are pushed to larger apocenter orbits. This observation invites us 
to consider a toy model for SIDM halos where the effect of SIDM 
is confined to a region (smaller than a radius of about rb) wherein 
particles redistribute energy and move towards a constant density 
isothermal core. We will develop this model further to explain the 
scaling relations between core size and halo mass in jj6] 

The circular velocity curves for the same set of halos discussed 
above are shown in Figure [5] The SIDM rotation curves rise more 
steeply and have a lower normalization than for CDM within the 
NFW scale radius r s . This brings to mind the rotation curves ob- 
served for low surface brightness galaxies and we will explore this 
connection later. Note though that the peak circular velocity Vmax 
actually is slightly higher for the SIDMi case because of the mass 
rearrangement (evident in the density profiles in Figure |4j briefly 
discussed in the last paragraph. At radii well outside the core ra- 
dius, the rotation curves of the CDM and SIDMi halos converge, 
though this convergence occurs beyond the plot axes > r s for most 
of the halos shown. 

An appreciation of why the density profiles of SIDM halos 
become cored can be gained from studying their velocity disper- 
sion profiles compared to their CDM counterparts, as illustrated 



in Figure [6] Here v ims is defined as the root-mean-square speed 
of all particles within radius r. While the CDM halos (black) are 
colder in the center than in their outer parts (reflecting a cuspy 
density profile) the SIDM halos have hotter cores, indicative of 
heat transport from the outside in. Moreover, the SIDM halos are 
slightly colder at large radii, again reflecting a redistribution of en- 
ergy. As discussed in the introduction, it is this heat transport that 
is the key to understanding why CDM halos differ from SIDM ha- 
los in their density structure dBalberg et alj2 002; Coli n et alj2002t 
I Aim & Shapiro! 120051 ; iKoda & Shapiro! 1201 ll) . The added thermal 
pressure at small radii is what gives rise to the core. The SIDMi 
simulations have sufficient interactions that they have been driven 
to isothermal profiles for r/r s < 1, while for SIDMo.i the v ims 
profiles typically begin to deviate from the CDM lines only at 
smaller radii, r/r s ~ 0.2, reflecting the relatively lower scatter- 
ing rate. 

The deviations in the SIDM ti rms profiles compared to CDM 
appear to set in at approximately the radius where we expect every 
particle to have interacted once in a Hubble time. This is explored 
directly in Figure[7j where we present a proxy for the local scatter- 
ing rate as a function of distance from the halo center: 



p(r)Vrms(r) oc T(r) (cr/mY 



(13) 



We have divided out the cross section so it is easier to compare 
the SIDMo.i and SIDMi cases. FigureUJpresents this rate proxy in 
units of 1 Gyr cm 2 /g: for the SIDMi case (with a/m — 1 cm 2 /g) 
the radius where a typical particle will have scattered once over a 
10 Gyr halo lifetime is p(r)w rms (r) = 0.1. For the SIDMo.i case 
(with a/m — 0.1 cm 2 /g), the ordinate needs to be ten times higher 
(~ 1) in order to achieve the same scattering rate. 

By comparing Figure [7j to Figure [6] (and to some extent to 
all Figures 4-6) we see that the effects of self-interactions do be- 
come evident at radii corresponding to pv ims ~ 0.1 for SIDMi (at 
r/r B ~ 0.8) and pv vmE ~ 1 for SIDMo.i (at r/r E ~ 0.2). Inter- 
estingly, for the SIDMi halos this interaction radius is fairly close 
to the Burkert scale radius (shown by the blue arrows). It should be 
kept in mind, however, that the structure of halos can be affected 
to larger radii because particles scattering in the inner regions can 
gain energy and move to larger orbits. A careful inspection of the 
density and rotation velocity profiles shows that this is indeed the 
case. 

We will discuss these findings in more detail in Sections [6] 
and UJ In particular, in ij7] we present an analytic model aimed at 
understanding how the central densities and scale radii of SIDM 
halos are set in the context of energetics. But before moving on to 
those issues, we first explore halo substructure in SIDM. 

5.4 Substructure 

The question of halo substructure is an important one for SIDM. 
One of the original motivations for SIDM was to reduce the num- 
ber of subhalos in the Milky-Way halo in order to match the relativ e 
dearth of observed satellite galaxies (Spergel & Stei nhardHl2000l) . 
However, the over-reduction of halo substructure is now recog- 
nized as a negative feature of SIDM compared to CDM, given the 
clear evidence for galaxy- size subhalos throughout galaxy clus- 
ters dNataraian et alj|2009h and the new discoveries of u ltra-faint 
galaxi es around the Milky Way (see Willma^ d2010h and lBullock] 
(2010) for reviews). In fact, one of the most stringent constraints 
on the self-interaction cross section comes from analytic subhalo- 
evaporation arguments dGnedin & O striker 200l|). 

Figure [8] demonstrates that the effects of subhalo evaporation 




Figure 4. Density profiles for our six example halos from our SIDMi (blue stars) and SIDMq.i (green triangles) simulations and their CDM counterparts. 
With self-interactions turned on, halo central densities decrease, forming cored density profiles. Solid lines are for the best NFW (black) and Burkert (blue) fits, 
with the points representing the density at each radial bin found by AHF. The arrow indicates the location of the Burkert core radius r^. r B is the NFW scale 
radius of the corresponding CDM halo density profile (black solid line). Burkert profiles provide a reasonable fit to our SIDMi halos only because j*b ~ r s 
for a/m = 1 cm 2 /g, so a cored profile with a single scale radius works. As discussed in [Qthis is not the case for a/m = 0.1 cm 2 /g and thus Burkert 
profiles are not a good fit to our SIDMq.i halos. 



in SIDM are not as strong as previously suggested on analytic 
grounds. Here we show the cumulative number of subhalos larger 
than a given Knax for a sample of well-resolved halos in our CDM 
(solid), SIDMo.i (dotted), and and SIDMi (dashed) simulations. 
The associated virial masses for each host halo are shown in the 
legend. The left panel presents the Vmax function for all subhalos 
within the virial radius of each host and the right panel restricts the 
analysis to subhalos within half of the virial radius. We see that gen- 
erally the reduction in substructure counts at a fixed V mxc is small 
but non-zero and that the effects appear to be stronger at small radii 
than large. Similarly, there appears to be slightly more reduction of 
substructure in the SIDM cluster halos compared to the galaxy size 
systems. 

We can understand both trends, 1) the increase in the differ- 
ence between the CDM and SIDM Vmax functions as Afvtr in- 
creases and 2) the increase in the difference as one looks at the 
central regions of the halo, using the results from the previous sec- 
tion as a guide. The typical probability that particle in an SIDM 
subhalo will interact with a particle in the background halo is 

P « (Phost{i-){v/rn)v orb (r)) T T, (14) 

where v or b(r) is the orbital speed of the subhalo at position r, phost 
is the mass density of the host halo, and T is the orbital period. 
The typical speed of the subhalo is similar to the rms speed of 
the smooth component of the halo, and thus phost(j){a /m)v rb{^) 



should be similar to the function we show in Figure [7] At fixed 
r/r B we expect P to scale with V max as Vmax/ r max (given that 
p s oc V^ a x/r^ a x), which is a very mildly increasing function 
of Vmax over the range of halo masses we have simulated. Note 
though that we expect scatter at fi xed halo mass beca use of the 
scatter in the Vmax — r max relation (Bulloc k~et alj200l1) . 

While the increase in destruction of subhalos with host halo 
mass is not strong, it is clear from the above arguments that subha- 
los in the inner parts of the halo (r/r s <C 1) should be destroyed 
but the bulk of the subhalos around r/r s ~ 1 and beyond should 
survive for a/m = lcm 2 /g. This effect is strengthened by the 
fact that subhalos in the innermost region of the halo were accreted 
much longer ago than subha los in the outskirts , so they have ex- 
perienced many more orbits jRocha etalj|201ll) . These arguments 
explain the comparisons between the subhalo mass functions plot- 
ted in Figure[8] Our arguments demonstrate that a large fraction of 
the subhalos found in CDM halos (most of which are in the outer 
parts) would still survive in SIDM halos for a/m values around or 
below 1 cm 2 /g. 

Overall in the previous two sections we have seen that the effects 
of self-interactions between dark matter particles in cosmological 
simulations are primarily in the central regions of dark matter ha- 
los, leaving the large scale structure identical to our non-interacting 
CDM simulations. Thus we retain the desirable features of CDM 
on large scales while revealing different phenomenology near halo 




Figure 5. Circular velocity profiles for our example selection of six well resolved halos from our CDM, SIDMi and SIDMo.i simulations. The magnitude of 
the circular velocity at small radii r < r s is lowered for all halos when self-interactions are turned on. r s is the NFW scale radius of the corresponding CDM 
halo density profile. 



centers. In the following section we will move to explore how the 
properties of SIDM halos presented here scale with halo mass. 



6 SCALING RELATIONS 

In the previous section we saw that while SIDM preserves the CDM 
large scale properties of dark matter halos, self-interactions in the 
central regions of halos result in a decrease of central densities and 
the formation of cores in their density profiles. We found that the 
density profiles of halos from our SIDMi simulations can be rela- 
tively well fit by Burkert density profiles inside r ~ 2 — 3r s (see 
Figure 0- Here we define a sample of well resolved halos from all 
our SIDMi simulations and use Burkert fits to their density pro- 
files in order to quantify their central densities and core sizes. We 
then provide scaling relations of dark matter halo properties with 
maximum circular velocity V max - 

The sample of halos used for the rest of this section consists 
of the two host halos in our SIDMi-Zll and SIDMi -Z12 simula- 
tions together with the 25 most massive halos from our SIDMi - 
50 and the 25 most massive halos from our SIDMi-25 simula- 
tions. That gives us a total of 52 halos spanning a range Vm ax = 
30 - 860 km/s or A/ vir = 5 x 10 11 - 2 x 10 14 M . For this set 
of halos the innerm ost resolved radius, defined by Equation 20 in 
IPower et all J2003h . is always smaller than one third of the Burkert 
scale radius from which we define the si zes of cores. It is v ital that 
we do a conservative comparison to the lPower et alj {2003) radius 
because both gravitational scattering and self-interactions lead to 



the same phenomelogical result of constant density cores. Most of 
the halos (other than the 52 we select here) do not pass this test 
well enough for the core set by self-interactions to be resolved with 
confidence. This desire to be conservative in our presentation of 
scaling relations forces us to find these relations from only a small 
sample of halos for SIDMi and leaves us with basically no halos 
to find scalings for SIDMo.i. Also one has to keep in mind that 
our SIDMi relations could be biased by selecting only the most 
massive halos in our full box simulations. Evidently higher res- 
olution simulations are necessary to find definitive answers. It is 
reassuring however that the scaling relations derived from our ana- 
lytical arguments in ij7]agree so well with the ones presented here 
for a/m — 0.1 cm 2 /g. 

We have checked that for all of our halos we resolve the scat- 
tering rate out to at least four times the Burkert scale radius. Outside 
of this point the scattering rate is underestimated because of our 
choice of the self-interaction smoothing length relative to the inter- 
particle spacing (see iJ3). However, the expected scattering rate is 
negligible with respect to the Hubble rate outside that radius (Fig- 
ure[7). Moreover, we have re-run our 50/i _1 Mpc boxes for a range 
of SIDM smoothing values and found identical results. Thus we 
consider our sample to be well resolved. 

Eight halos in our sample are undergoing significant interac- 
tions and have density profiles that are clearly perturbed even in the 
CDM runs. We include these eight systems in all of the following 
plots but indicate them with open symbols. We do not use them in 
the best fits for the scaling relations that we provide. 

We start by examining the global structure of halos as char- 
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Figure 6. Velocity dispersion profiles for our six example halos from our SIDMi and SIDMo.i simulations over-plotted with their CDM counterparts. The 
velocity dispersion is inflated at small radii and slightly suppressed at large radii. The effects set in at approximately the radius where SIDM particles experience 
at least one interaction on average over the lifetime of the halo (see Figure|7J. 



acterized by the maximum circular velocity V ma x and the radius 
where the rotation curve peaks, r max . The relationship between 
Vmax and r max provides a simple, intermediate-scale measure of 
halo concentration and we aim to investigate any differences be- 
tween SIDM and CDM. Figure [9] shows the Vm ax — r max rela- 
tion for CDM (black) and SIDMi (blue) halos. We can see that 
small differences of about 10% exists in both V ma x and r max , with 
SIDMi halos having larger values for Vm ax and smaller for r max . 
This was already evident in Figure [5] where the circular veloc- 
ity curves of SIDMi halos seem to peak at slightly smaller radii 
and slightly larger velocities than their CDM analogs, even though 
SIDMi curves decrease more steeply at the center. 

The apparent difference is consistent with a picture where en- 
ergy exchange due to scattering redistributes the SIDM dark matter 
particles, with many of the tightly bound particles scattered onto 
less bound, high apocenter orbits. Since the radius at which self- 
interactions are significant (see Figure[7} is smaller than (but close 
to) r s , it is entirely reasonable that the scattered particles lead to 
a new r max for SIDMi that is smaller than the CDM r max and a 
Vm ax that is larger. Notice that the slope of the Vmax— '"max relation 
is unchanged from CDM to SIDMi . The best-fit relations are: 
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(15) 



We continue this discussion by considering the sizes of cores 
in our SIDMi simulations as a function of Vm ax - The core sizes of 



halos are quantified by the scale radius in the Burkert fit to their 
density profiles, namely r b in Equation [T2] Figure [Tol shows that 
for this relation a single power law holds along the whole range 
of our sample. We will come back to this result in our discussion 
section (ij8} on extrapolating to smaller and larger V max values to 
make contact with observations of cores in galaxies and clusters. 
The power law that best fits our data is given by 
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If we fit to M v i r instead of V ma x we get 



r b = 2.21 kpc 
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(16) 



(17) 



We note that the scaling with V max is close to that expected 
for r max or r s . We show this explicitly by fitting for the core size 
of SIDMi halos r b as a function of the NFW scale radius r s of 
their CDM counterparts, as shown in Figure [T71 We find that the 
ratio of the core size of a SIDMi halo to the scale radius of the 
corresponding CDM halo varies very mildly with Vm ax - In other 
words, the core sizes are a fixed fraction of the CDM halo scale 
radius. The relation that best fits our data is given by 



rb 



0.71 



10 kpc 



(18) 



This underscores the point that r b and r s are closely tied to each 
other and the fact that they are numerically so close to each other 




Figure 7. Estimate of the local scattering rate modulo the cross section pv TulB = r((x/m) — 1 for six well resolved halos from our CDM, SIDMo.i, and 
SIDMi simulations. The quantity is scaled by 1 Gyr cm 2 /g, such that 1 in these units means that each particle has roughly one interaction per Gyr in SIDMi 
and 0.1 per Gyr in SIDMo.i . Based on this argument, the effects of self-interactions in the properties of halos over ~ 10 Gyr should start to become important 
when the ordinate is greater than about 0.1 in SIDMi (r/r B ~ 0.8) and greater than about 1 in SIDMo.i (r/r B ~ 0.2). Comparisons to Figures 4-6 indicate 
that this is indeed the case. 



is the reason why a cored profile with a single scale (like a Burk- 
ert profile) provides a reasonable fit to our SIDMi halos. We will 
explain this striking behavior using an analytic model in the next 
section. 

The central densities in SIDMi halos can be defined either as 
the Burkert profiles scale density or as the density at the innermost 
resolved radius. We have found that both definitions give similar 
results with no significant differences. In Figure[l2] we show how 
the Burkert scale density pb scales with V max . The trend in the 
Pb — Vmax relation is not as strong as for the — Vmax relation, 
with a scatter as large as about a factor of 3. We will come back 
to the implications of this result in our discussion section (fj8). The 
relation that best fits our data is given by 

Pb = 0.015 Me/PC? (-^J (19) 

If we fit to M v i r instead of Vmax we get 

/ m- \-° 19 
p b = 0.029 M /pc 3 (toI%) • (20) 

We urge caution when using the above fits to the central densities 
as it is likely to be affected by our small sample size given the 
large scatter. The toy model discussed in the next section predicts 
a slightly stronger scaling with Vmax • However, the typical densi- 
ties of order 0.01 M /pc 3 for galaxy halos and 0.001 M Q /pc 3 for 
cluster halos (see Figure [T2l are in line with the predictions of the 
analytic model. 



In this section we have presented scaling relations for the properties 
of halos in our SIDMi simulations. Our limited resolution allows 
us to use only 52 halos spanning a modest mass range, from which 
we throw out eight systems that are undergoing mergers. Admit- 
tedly, this sample is not large enough to be definitive, especially 
in regards to scatter. However, the strong correlation between the 
SIDM core radius r b and the counterpart CDM scale radius r s is 
clearly statistically significant and the general trends provide a use- 
ful guide for tentative observational comparisons - a subject we 
will return to in the final section below. 



7 ANALYTIC MODEL TO EXPLAIN THE SCALING 
RELATIONS 

In this section we develop a simple model to understand the scal- 
ing relations shown in Sj6] This model is based on identifying an 
appropriate radius ri within which self-interactions are effective 
and demanding that the mass as well as the average velocity dis- 
persion within this radius is set by the mass and the average ve- 
locity dispersion (within the same radius) of the same halo in the 
absence of self-scatterings. The mass loss due to scatterings in the 
core should be insignificant because particles rarely get enough en- 
ergy to escape and this implies that the mass within n should be 
close to what it would have been in the absence of self-interactions. 
This also implies that the potential outside n is unchanged from its 
CDM model prediction, but tends to a constant value faster inside 
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Figure 8. Subhalo cumulative number as a function of halo peak circular velocity (V ma x) for several well-resolved halos in our CDM (solid), SIDMo.i 
(dotted), and SIDMi (dashed) simulations. When looking at all subhalos within r < r v ; r (left), the differences are small and the slope of the subhalo V m ax 
function is the same for the CDM and SIDM cases. The offset in the subhalo Vmax function increases when we look only at subhalos inside r < 0.5 r v ir 
(right panel), showing that SIDM suppresses the number of subhalos in the central regions of halos more strongly. 



n. Within this set of approximations, the dominant effect due to 
scatterings is to re-distribute kinetic energy in the core, while keep- 
ing the total kinetic energy within n the same as it would have had 
before self-interactions became important. We have looked at the 
kinetic energy profiles in the best-resolved halos in our simulations 
and have confirmed that this is indeed a good approximation. Note 
that in this picture, there is a clear demarcation of time-scales such 
that the inner halo structure (say r < r s ) is set (the same way as in 
CDM model) well before self-interactions become important. For 
cross sections much larger than what we are interested in here, this 
need not hold. 

To set up the model, we start by recalling that self-interactions 
work to create an isothermal core (see Figure [6} that is isotropic 
(both spatially and in velocity space). Using the spherical Jeans 
equation, one can then see that for a system with these properties 

2nC 1 Gp(0)r 2 o , 



2 



(21) 



where we have defined ro to be the expansion parameter such that 
p{r)a r (r) 2 = p(0)oy(0) 2 (l - £(r/r ) 2 ) when r < r , and a r 
is the radial velocity dispersion. The form of the Taylor expansion 
for p(r)a r (r) 2 is dictated by the Jeans equation for density pro- 
files that tend to a constant value, as may be readily ascertained 
by taking the derivative of p(r)a r (r) 2 . To fix ro, we will choose 
it to be equivalent to the Burkert scale radius where the density is 
one-fourth of the central density. The parameter £ encapsulates un- 
certainties from the profile and velocity dispersion anisotropy in the 
outer parts of the halo. We test various models and find that a range 
of 2-3 for £ is largely consistent with most parameterizations and 
hence we fix £ = 2.5. If we specify the central velocity dispersion, 
then with an additional constraint on the core region (i.e., ri), we 
would be able to back out both the core radius and the core density. 
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Figure 9. r m ax vs. Vmax for our combined sample of well resolved halos 
from our SIDMi and CDM simulations. Open symbols correspond to halos 
for which the density profiles showed signs of being perturbed, thus they 
were not included in the best fit of the relation. Small differences of about 
10% exists in both Vmax and r m ax, however the slope of Vm 
lation is unchanged from CDM to SIDMi. 



We then set v rm s,o equal to the average velocity dispersion 
squared (i.e., two times kinetic energy divided by mass) within 
the region n in the absence of self-interactions. This basically de- 
mands that the kinetic energy within r\ is unchanged from the value 
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Figure 10. Burkert scale radius vs. V max for our combined sample of 
well resolved halos from our SIDMi-50 (blue circles), SIDMi-25 (green 
stars), SIDMi-Z12 (cyan square) and SIDMi-Zll (red triangle) simu- 
lations. Open symbols correspond to halos that are undergoing mergers. 
These perturbed halos were not included in the fit for the scaling relation. 
A single power law holds along the whole range of our sample, suggesting 
that this dependence continues towards smaller and larger Vm ax values. 



it would have had in the absence of self-interactions. Note, how- 
ever, that we are setting the average velocity dispersion squared 
equal to v 2 ms> o and not the corresponding average in the SIDM 
halo. This is an approximation, but one that is degenerate with 
choosing the £ parameter. 

To finish specifying this model, we need a density profile for 
the region inside r 1 . A Burkert profile has a velocity dispersion pro- 
file (assuming isotropy) that asymptotes very slowly to the central 
dispersion. For small radii, the radial dispersion profile is slowly 
increasing (with radius) because of the r/Yb term in the Taylor 
expansion for the density profile. If we want a flatter central dis- 
persion profile (as is observed for the SIDMi halos), we can fix 
this by either assuming an isothermal profile or something like 
1/(1 + 1.52(r/r ) 2 ) 3/2 . The final results turn out to be qualita- 
tively similar for these profiles. Hence we adopt a Burkert profile 
for ease of comparison to the fits presented here and then check the 
results with more appropriate profiles later. Our two constraints (on 
the radial velocity dispersion and mass) fully specify the density 
and radial scales of the Burkert profile. 

In order to obtain scaling relations we need to estimate n, 
which demarcates the inner region where self-interactions are ef- 
fective from the outer region that is mostly undisturbed by the self- 
interactions. In reality, this divide will not be sharp but we will see 
that the main features of the scaling relations are well-captured by 
this simple model. We define ri to be the region where each par- 
ticle on average suffers one interaction. Since the region outside is 
assumed to be unperturbed by interactions, we may estimate n as: 

r(ri)t ag G = l-3/9CDM(ri)ui m8 ,cDM('*i) — ia gc = 1 , (22) 

m 

where we set age (i age ) to be 10 Gyr for now, keeping in mind 
that larger halos have a shorter age and that major mergers can re- 
set the timer. We will consider what happens when i agc is a func- 
tion of halo mass shortly. The factor 1.3 is (\v — u\)/ \f (v 2 ) for 
a Maxwellian distribution where u and v are the velocities of the 



two interacting dark matter particles. We have not attempted to use 
a more realistic velocity distribution since the dependence of this 
factor on a possible high-velocity cut-off to the distribution func- 
tion was found to be fairly mild. 

For the density profile in the absence of self-interactions, 
we assume a NFW profile and to fix the velocity dispersion 
we use the observed fact that the phase sp ace density is a 
power-law in radius dTavlor & Navarro] l200ll) . By noting that 
Urms,CDM(r) = (pcdm (r)/Q(r ) ) 1/3 and using a phase-space 
density profile Q ( r) = Q(r s )(r /r f )- v dTavlor & Navarroll200ll; 



Rasia ct alj |2004 lAscasibar etaU |2004 iDetmen & M cLaughlin 



20051 : lAscasibar & Gottl6berll2008l) . we may fully specify the de- 
pendence of n on the cross-section and halo parameters (say V max 
and r max ). For the phase-space density profile we use a power-law 
index 77 = 2 and Q(r B ) — 0.3/(GT^ nax r 2 nax ) derived from jointly 
fitting our relaxed CDM halos; these parameters are very similar to 
the fits provided in lAscasibar & G ottlober (2008). 

Let us first look at how n scales with r s in the NFW den- 
sity profile. One notes that p B = 1.72V^ ax /(Gr 2 nax ) and hence 
p s Kn ax oc V I 3 lax /?"m ax which is a very mildly increasing func- 
tion of Vmax as our Equation [T5] shows. Thus Equation [22] im- 
plies that ri/r s should be roughly a constant. Numerically, we 
find that ri/r 8 — 0.7 — 0.8 over the range of Vm ax of interest 
for a/m — 1 cm 2 /g. 

Having now specified ri, we are ready to look at the scalings 
of rh and pb. For our assumed value of f, v 2 ms ,o — 2.5Gpbr 2 . 
Thus we are looking for the value of rb /r s that solves, 



(r/r,)" 



(r/r a )(l + r/r s ) 2 Q{r s ) 



(ri) = (^rms.o)" 



(23) 



with the constraint that Mbiri) = AfNFw(ri) where Mnfw(^) 
and Mh(r) are the masses enclosed within radius r for NFW and 
Burkert profiles, respectively. We note that if rb/r s is not a strong 
function of Vmax and since we know ri/r s is a mild function of 
V max , then the mass constraint essentially sets phr^ / {p s r'i) to be 
a constant or p b r 3 oc (p s r 3 ). This implies « 2 ms , rb oc F max r max . 
Now Equation|23]sets v rms ,o ~ V m &x because n /r s is a mild func- 
tion of V max and it therefore follows that rb oc r s is a consistent 
solution to the above equations. As a check we note that assuming 
ri/r s = 0.7 — 0.8 gives » lmS) o ~ l.lVm ax , in reasonable agree- 
ment with our SIDMi simulation results (see Figure[6j. This sim- 
ple model thus predicts that rb /r s should not vary much with V max 
in agreement with the observed scaling relations from the SIDMi 
simulation. 

In detail, the model predicts that rb/r s = 0.5 — 0.6 for dwarf 
to cluster halos in good agreement with the fits to our SIDMi 
halos, but about 25% smaller for Vm ax ~ 100 km/s. It departs 
from the results of the simulation in predicting that rb/r s in- 
creases gently with V max , whereas Figure QT] predicts that this ra- 
tio should decrease gently with Vm ax . We find that this departure 
from simulations is likely related to the assumption of a constant 
age for all halos. To generalize our model, we use the results of 
IWechsler et alj «|2002) who show that the virial concentrations of 
halos are correlated with their formation times, and in particular 
Cvir = 4.1(1 + Zform) for a particular definition of formation time. 
We invert this equation to derive an estimate of the halo age us- 
ing Zf orm . With the age thus specified in Equation [22] we find that 
now rb/r s decreases gently with V max in substantial agreement 
with the fit to our simulations. Thus the reason that larger halos 
have a smaller rb /r s is because self-interactions have had less time 
to operate. We note that the values for the core radius in the an- 
alytic model with halo mass dependent t ago are uniformly about 
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Figure 11. Burkert scale radius in SIDMi halos vs. the NFW scale radius 
in their CDM counterparts. Points and labels are the same as in Figure [Tol 
There is a one-to-one correlation indicating that the core size of SIDMi 
halos scales the same as the scale radius of CDM halos with Vmax 

25% smaller, but this should not be a cause for concern given the 
approximation in demanding a sharp transition at r 1 . 

Given the Burkert core radius r b and the central velocity 
dispersion « r ms,o, one can easily check that the central density 
Pb is about O.OIMq/pc 3 for V ma x = 300km/s halos and 
0.005 M Q /pc 3 for Vmax = 1000 km/s in this analytic model. 
These numbers and the scaling with Vmax for pb (when including 
the halo mass dependent t agc ) are in good agreement with the den- 
sities in Figure[T2]and the fit in Equation[T9] As we have indicated 
before, the scaling relation for the central density should be inter- 
preted with care given the large scatter. Given the tight correlation 
between core radius and r s , it is possible that the substantial scatter 
in the central density arises in large part due to the scatter intro- 
duced by the assembly history in the concentration-mass relation. 
This has important implications for fitting to the rotation veloc- 
ity pr ofiles of low-surface brightness spirals dKuzio de Narav et al.l 
2010) and deserves more work. 

The simple model constructed above also provides insight into 
the core collapse time scales. In particular, as long as the outer part 
(region outside ri) dominates the potential well and sets the aver- 
age central temperature (or the total kinetic energy in the core), we 
do not expect core collapse. This is simply because core collapse 
requires uncontrolled decrease in temperature, which is prohibited 
here. Once r± moves out well beyond r max or to the virial radius, 
there is significant loss of particles and core collapse may occur if 
there are no further major mergers. The time scale for this process 
is much longer than the age of the universe for a/m = 1 cm 2 /g 
because the inner core is at n < r s after 10 Gyr for this self- 
interaction strength and we see no evidence for significant mass 
loss. 



8 OBSERVATIONAL COMPARISONS 

The goal of this section is to discuss our results in comparison to 
observationally inferred properties of dark-matter density profiles. 
In particular, we will focus on the core densities and core sizes. £13- 1 1 



Figure 12. Burkert scale density vs. Vmax. Points and labels are the same 
as in Figur JlOl The trend in the pb — Vmax relation is not as clear as for 
the r'b — Vmax relation, with a scatter of up to a factor of 3. 



presents our expectations for SIDMi and SIDMo.i . Our predictions 
for a/m — 1 cm 2 /g are anchored robustly to our simulations, 
though they do require some extrapolation beyond the mass range 
directly probed by our simulations (Vmax = 130 — 860 km/s). 
For a/m = 0.1 cm 2 /g the predictions are much less secure be- 
cause the associated core sizes are of order our resolution limit, 
thus we rely on our our analytic model more directly here. In i]8.2l 
we discuss our predictions in light of observations of dark-matter 
halos for a wide range of halo masses. In £18.31 we discuss our re- 
sults on subhalos in the context of past work and constraints on 
SIDM based on subhalo properties. 

Before proceeding with this discussion we would like to clar- 
ify how we quantify core sizes. In this work, we have fit the 
a/m = lcm 2 /g halos with Burkert density profiles. However, 
many observational constraints on cores on galaxy scales come 
from fittin g pseudo-isothermal density profiles with core size r p i to 
data (e.g. JSimon e t al. 2005: IKuzio de Narav et al . 2008), al though 
some constraints do come from Burkert modeling JSalucci et all 
l2012h . We found that pseudo-isothermal density profiles also give 
good fits to the inner regions of the SIDMi halos, but Burkert fits 
are better because of that profile's p oc r~ 3 dependence at large 
radii. For a pseudo-isothermal density profile (oc l/(r 2 + r 2 )), the 
density decreases to one-fourth the central density at 1.73 times its 
core radius r c . Thus, as a crude approximation, one may convert 
the Burkert radius to the equivalent pseudo-isothermal core radius 
by multiplying by a factor of 0.58 (r c ~ rt/1.73). 



8.1 Predicted Core Sizes and Central Densities in SIDM 

8.1.1 SIDM with a/m = 1 cm 2 /g. 

The central properties of dark-matter halos have been inferred from 
observations from tiny Milky Way dwarf spheroidal (dSph) galax- 
ies (Vmax ;$ 50 km/s) to galaxy clusters (Vmax > 1000 km/s). If 
we extrapolate the results from our set of SIDMi simulations using 
Eqs. l[l6j-{|20]l we predict that SIDM halos with a/m = 1 cm 2 /g 
would have the following (Burkert) core sizes and central densities: 
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For galaxy clusters (Vmax ^ 700 - 1000 km/s): 

r b ~ (95 - 155) kpc ; p b ~ (0.005 - 0.004) M pc~ 3 
For low-mass spirals (V ma x — 50 — 130 km/s): 

r b ~ (3 - 10) kpc ; p b ~ (0.02 - 0.01) M Q pc" 3 
For dwarf spheroidals galaxies (Vmax — 20 — 50 km/s): 

r b ~ (0.9 - 3) kpc ; p b ~ (0.04 - 0.02) M pc" 3 

Although we can not completely determine the scatter in our scal- 
ing relations due to low number statistics, it is important to note 
from Figs. [10] and [T2] that a scatter of at least a factor of 2 in core 
sizes, and at least a factor of 3 in central densities, is expected for 
a given V max . We suspect that these differences are in large part 
a result of the diversity of merger histories of dark-matter halos. 
Note that the V m ax-r m ax and Q(r s ) scalings assumed in the ana- 
lytic model are the median values. The strong dependence of the 
SIDM halo profiles on these quantities makes it clear that the scat- 
ter in these relations will introduce significant scatter in the halo 
core sizes and core densities. Thus the analytic model should also 
provide a simple way to understand (some of the) scatter seen for 
SIDMi halo properties. In future work we will characterize the re- 
lation between the core properties and merger history in the context 
of a detailed discussion of the scatter in the scaling relations, espe- 
cially on scales that we do not resolve with our current simulations. 

S. 1.2 SIDM with a/m = 0.1 cm 2 /g 

As discussed in i]5.3| our SIDM .i simulations are not well enough 
resolved to definitively measure a core radius for any of our ha- 
los, much less define scatter in that quantity. Nevertheless, our best 
resolved systems do demonstrate some clear deviations from CDM 
and allow us to cautiously estimate individual core densities. Refer- 
ring back to Figure 4, we see that in our two best resolved cluster 
halos (at M vlr 10 14 M ) the SIDM . i core densities approach 
~ 0.01 M /pc 3 - each at least a factor of ~ 3 denser than their 
SIDMi counterparts. Similarly, in our Z12 Milky- Way case, the 
SIDMo.i core density appears to be approaching ~ 0.1 M /pc 3 
compared to ~ 0.02 M /pc 3 in the SIDMi case. 

Given the lack of well-resolved halo profiles, it is worth ap- 
pealing to the analytic model presented in 5j7]to estimate core radii 
for SIDMo.i. Using exactly the same arguments (including the halo 
mass dependent age), we find that ri/r s ~ 0.05 — 0.12 in the 
100-1000 km/s Vmax range and a corresponding Burkert core ra- 
dius r%/r B ~ 0.09 — 0.17. We note that the Burkert radius is close 
to but slightly larger than n. It is important to keep in mind that 
in this analytic model we are only explicitly fitting the inner "self- 
interaction zone" of r < r\. This does not imply that the entire halo 
has to be well-fit by the Burkert profile. Recall that a single-scale 
Burkert profile only works as well as it does for a/m = 1 cm 2 /g 
because rb ~ r 3 , such that to a good approximation there is only 
one relevant length scale. For the smaller cross section that we 
are now considering we expect the core and NFW scale radii to 
be widely separated, suggesting that a generic functional form for 
SIDM halos should have two scale radii. A wide separation be- 
tween the SIDMo.i core and r s does appear to be consistent with 
the highest resolved halos presented in Figure 4. However, we note 
that given the strong correlation between r b /r B , we still expect a 
one-parameter family of models for a given a/m. 

To see how dependent our results are on the shape of the in- 
ner halo profile, we modify the analytic model to include a density 



profile that decreases with radius as 1/(1 + (r /r c ) 2 ) a ' 2 . For this 
density profile, the velocity dispersion profile has the right form to 
match our simulation results. The price we pay is the introduction 
of a new parameter a. We set this parameter alpha by additionally 
demanding that the slope of the mass profile (i.e., density) is con- 
tinuous at n , so that the mass profile joins smoothly with the NFW 
mass profile. This picks out a narrow range a — 5.5 — 7.0 as the 
solution over most of the Vmax range of interest (with smaller val- 
ues corresponding to lower V ma x). Interestingly, this implies that at 
n, the slope of the density profile is very close to —2 for the entire 
range of Vmax values of interest. Note that while the mass profile is 
continuous, the slope of the density profile is not matched smoothly 
at n (since the slope of the NFW profile would be closer to —1 at 
fi <C r B ) . This probably signals that if the matching were not done 
sharply (at ri), the density profile of SIDM would overshoot that 
of CDM and catch up at some radius beyond n (as is seen in the 
comparison of SIDMi and CDM density profiles). 

As a check we apply this a-model to a/m — 1 cm 2 /g case 
and find that the results are qualitatively the same as the model 
with the Burkert profile. The quantitative differences are at 20% 
level with the densities being smaller and inferred Burkert core 
radii (where density is 1/4 of the central density) larger compared to 
the Burkert profile model. The predicted slope of the density profile 
at n is close to —2.5 implying a smoother transition to the NFW 
profile (since ri ~ r B for a/m — 1 cm 2 /g), as is seen Figure|4] 

For the a/m — 0.1 cm 2 /g case, we obtain r c /r B = 0.08 — 
0.17 and an equivalent Burkert core radius (where the density is 
one-fourth of the central density) r b /r 3 — 0.06 — 0.14, in sub- 
stantial agreement with the results we obtained using the Burkert 
profile. Thus our analysis would suggest that core sizes ~ 0.1r s 
for a/m = 0.1 cm 2 /g. The results from the analytic model for 
a/m — 0.1 cm 2 /g also seem consistent with our simulation re- 
sults; see Figure 6 where the u rm s profiles for SIDMo.i start to 
deviate from CDM at ~ 0.2r s . 

Based on the discussion above we conclude that for a/m = 
0.1 cm 2 /g we expect: 

For galaxy clusters (Vmax - 700 - 1000 km/s): 

r b ~ (16 - 20) kpc ; p b ~ 0.04 M pc" 3 
For low-mass spirals (Vmax — 50 — 130 km/s): 

r b ~ (0.6 - 2.5) kpc ; p b ~ 0.2 - 0.1 M pc -3 
For dwarf spheroidals galaxies (Vmax — 20 — 50 km/s): 

r b ~ (0.2 - 0.6) kpc ; p b ~ 0.5 - 0.2 M pc" 3 

These values do not include the scatter from mass assembly his- 
tory. It is probably reasonable to assume a factor of 2 scatter for 
both core radii and core densities based on what we see in SIDMi. 
It is also possible that the core densities are ~ 50% smaller than 
what we would see in simulations, given that the SIDMi simula- 
tions have core densities that are somewhat larger than the predic- 
tions from the analytic model. For the dwarf spheroidal galaxies, 
the values should be interpreted with caution as it is the prediction 
for field halos with Vmax range 20 — 50 km/s. 

While these values are somewhat tentative compared to those 
presented above for SIDMi (given our lack of direct simulation 
fits), two factors are reassuring. First, the analytic model is based 
on the simple assumption that scattering redistributes kinetic en- 
ergy within the inner halo and the non-trivial aspect of the model is 
defining this "inner halo" region. There is no reason to suspect that 
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this assumption or the prescription breaks down for SIDMo.i halos 
when it works so well in describing the SIDMi halos. The pre- 
dicted densities are in line with those inferred for the best resolved 
halos in our SIDMo.i simulations (shown in Figure 4 and discussed 
above). For the core radii, we reiterate that the label "rb" should be 
interpreted (according to its definition in the analytic model) as the 
radius where the density reaches one-fourth the asymptotic core 
density. The overall profile of a halo with such a small core com- 
pared to r s will not be fit by the Burkert form. Note that the strong 
correlations we predict between the core radius and the NFW scale 
radius raise the intriguing possibility that the SIDM halos may be 
also well fit (modulo scatter) by a single parameter profile as is the 
case for CDM. 

Next, we compare our predictions for SIDM core properties 
against data and show that the core radii and densities appear to be 
consistent with that seen in real data, motivating future simulations 
with high enough resolution to resolve cores in SIDMo.i halos. 



8.2 Observed Core Sizes and Central Densities vs. SIDM 

In this section, we explore the predictions for the properties of den- 
sity profiles with SIDM in the context of observational constraints 
on density profiles. We also revisit previous constraints on SIDM 
from observations in light of our simulation suite. 



8.2.1 Clusters 

One of the tightest SIDM constraints from the first generation of 
SIDM studies emerged from one cluster simula t ion a nd one ob- 
served galaxy cluster. Specifically. I Yoshida et al.l ( feOOOh simulated 
an individual galaxy cluster with different SIDM cross sections. 
When comparing the core size of this simulated cluster to the 
core size estimated bv lTvson etal] Jl998h for CL 0024+1654, they 
found that the observed core in CL 0024+1654 would be consistent 
with SIDM only if cr/m < 0.1 cm 2 /g. Since that time, evidence 
has emerged that this particular cluster is undergoing a merger 
along the line of sight JCzoske et al.l200lll2002l : IZhang et al.l2005l ; 
IJee et alj|2007l : Heboid : lUmetsu et ai1l2010l) . Thus, this cluster is 
not the ideal candidate for SID M constraints base d on the proper- 
ties of relaxed halos, and the lYoshida et alj J200C|> constraint is not 
valid in this context. 

Using X-ray emission, weak lensing, strong lensing, stellar 
kinematics of the brightest cluster galaxy (BCG) or some combi- 
nation thereof, the mass distributions within a nu mber of galaxy 
cluster s have been mapped in the past decade. lArabadiis et al.l 
d2002l) placed a conservative upper limit of 75 kpc on the size of 
any constant-density core, and an average density within the inner 
50 kpc of ~ 0.025 Mq/pc 3 for an halo with an estimated mass 

M ~4x 10 14 M R . 

Sand et all b004h,ISand et alj J2008h .iNewman etal] J2009). 



and lNewman et al] d2011 ) all find central density profiles in clus- 



ters shallower than the NFW CDM prediction. The difference in 
the work between these authors and others is that they use stel- 
lar kinematics of the BCG to constrain the density profile of the 
cluster dark-matter halo on small scales. While this probe of the 
density profile is more sensitive on small scales than strong lensing 
is, proper inference of the dark halo properties depends on accu- 
rate modeling of the BGC density profile and equilibrium structure. 
They have typically assumed a "gNFW" profile in order to con- 
strain the central densities: p(r) oc l/(x 3 (l+x) 3 ~ 9 ) with r = xr B 
and the NFW form obtained when g — 1. The INewman et alj 



(2009) and lNewman et alj ( 1201 ll) mass models of M ~ 10 15 M Q 
clusters show average dark matter central densities within 10 kpc 
of ~ 0.03 - 0.06 M Q /pc 3 and r s of order 100 kpc. Note that 10 
kpc is typically the smallest radii our simulations can resolve. 

ISaha et alj d2006h and I S aha & Read J2009h studied the mass 
structure of 3 cluster halos from gravitational lensing and obtained 
density profiles that are consistent w ith p oc r -1 outsid e the 
inner 10 — 20 kpc r e gions . Similarly iMorandi et alj d2010l) and 
Mora ndi & Limousin] ( 20121) find that the radial mass distribution 
of cluster dark-matter halos are consistent with NFW predictions 
outside 30 kpc in projection. The CLASH multi-cycle treasury pro- 
gram on the Hubble Space Telecope is finding man y new strongly 
lensed galaxies in about a set of 25 massive clusters jPostman et alj 
l2012h . Initial results from this program show that the total density 
profile of these clusters (or total density minus the brightest cluster 
galaxy), if modeled as spherically symmetric, are consistent with 
NFW predictions for the halo alone if the gNFW functional form 
is use d in the fit Zitrin et alJl201ll;ICoe et al.ll2012l;IUmetsu et alj 
l2012h . However, IMorandi et all d2010h and lMorandi et all ( 1201 lh 
find that spherical mass modeling of galaxy clusters typically re- 
sults in an overestimate of the the cuspiness of the density profile, 
althou gh axially-symmetric modeling is found to lead to underesti- 
mates (Meneghetti et al. 2007). Thus, the present status of the den- 
sity profiles of the CLASH clusters is unclear and clearly an inter- 
esting data set to look forward to. 

We note here a complexity involved in using the lensing re- 
sults to constrain SIDM models. Lensing provides mass in cylin- 
ders along the line of sight and this 2D mass profile is sensitive 
to mass from a large range of radii. As an example, lets consider 
mass within 30 kpc in projection. If we were to do something ex- 
treme and create a zero density core inside 30 kpc sphere, the dif- 
ferences in the 2D mass profile would be less than a factor of 2 for 
clusters in the 10 14 " 15 M mass range. For SIDMo.i, the differ- 
ences are comparatively benign. Our analytic model predicts that 
differences relative to CDM at about 0.1r s (which is 10 — 40 kpc 
for 10 14 " 15 M virial mass range) are 20-30%, which implies 
SIDMo.i surface mass density profiles are very similar to CDM 
on these scales. But for SIDMi the expected differences would be 
measurably large. 

On a related technical note, we discourage the use of the 
gNFW functional form when thinking about models that deviate 
from the CDM paradigm. In the SIDM case, for a/m < 1 cm 2 /g, 
there will generically be two scale radii: one is the NFW-like 
scale radius which is the result of hierarchical structure formation 
(Lithwick & Dalalll20T ]]), and the second is the core radius from 
dark-matter self-scattering. For a/m = 1 cm 2 /g, as we explained 
in detail in jJJJ me two scales are about the same. If most of the clus- 
ter data constrain the density profile beyond a SIDM core, as they 
may for weak lensing and X-ray studies, the gNFW or NFW fit is 
dominated by those data, and a core will not be "detected" in the fit. 
In future work, we will simulate halos with a broader range of a/m 
and provide SIDM-inspired density profiles to the community. 

The results discussed above seem to suggest that the density 
profile beyond about 25 kpc should be close to the predictions from 
the NFW profile. To test this we plot the average physical den- 
sity within 25 kpc for well-resolved halos in our CDM (black), 
SIDMo.i (green), and SIDMi (blue) simulations in Figure[T3] We 
see that for the most massive halos, the a/m = 1 cm 2 /g run pro- 
duces densities at 25 kpc that are ~ 2 — 3 times lower than their 
CDM counterparts. Thus it seems like the measured densities in 
clusters rule out a/m — 1 cm 2 /g SIDM model. At the same time 
the cr/m = 0.1 cm 2 /g simulations are quite similar to CDM at 
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Figure 13. Dark matter average density within 25 kpc vs. V max for re- 
solved halos in our CDM, SIDMo.i, and SIDMi simulations. It is clear 
that SIDMi has significantly lower densities than CDM halos at group and 
cluster scales. For the SIDMo.i model, the differences are muted and only 
appear on cluster scales. Thus observations of central densities in clusters 
likely provide the most promising avenue to look for signatures of SIDM 
with cross sections in the vicinity of 0.1 cm 2 /g. 



these these radii, though beginning to show some differences as 
we discussed earlier in this section. Analyses that combine infor- 
mation from X-rays, lensing an d BCG stellar kinema tics seem to 
suggest lowered densities (e.g.. iNewman et all d201 if) ) that would 
be compatible with SIDMo.i. Given this outlook, it is reasonable to 
conclude that estimates of the central dark matter density in clusters 
will provide essential tests of interesting SIDM models. 



8.2.2 Low-Mass Spirals 

For low-mass spirals with maximum circular velocities in the 
range 50 — 130 km/s, constant-density cores with sizes of ~ 
0.5 — 8 kpc and central densiti es of approximately ~ 0.01 — 
0.5 Mp/pc 3 have been observed dde Blok et al J200ll : Simon et"al] 
20051: Isanchez-Salcedol l2005t iKuzio de Narav et aT]|2008l. |2oTcF 
Oh et alJl201ll : ISalucci et alj|2012h . Similar to what we found for 
clusters scales, SIDM with a/m = 1 cm 2 /g would be able 
to reproduce the largest core sizes observed in low-mass galax- 
ies but it predicts central densities that are too low. SIDM with 
a/m = 0.1 cm 2 /g would be much more consistent. Moreover, 
the predicted log-slope of the density profile at 500 parsecs for 
a/m — 0.1 cm 2 /g halos in the 50-130 km/s ra nge is —0.5 to , 
both facts consistent with results from THINGS dOh et alj|2oTTh . 
Note that the slope at 500 pc for the a/m = 1 cm 2 /g model is 
in the same V max range, which is not consistent with the scatter 
seen in the data. 

We conclude, as before, that the observed densities and core 
radii are not consistent with SIDMi but are fairly well reproduced 
in SIDM models with a/m ~ 0.1 cm 2 /g. 



8.2.3 Dwarf Spheroidals in the Milky Way Halo 

The least massive and most dark-matter-dominated galaxies pro- 
vide an excellent setting to confront the predictions of dif- 



ferent dark matter models w ith observations. Recent work by 
iBovlan-Kolchin et alj Jioi ld lbh has found that the estimated cen- 
tral densities of the bright Milky Way dwarf spheroidal satellites 
are lower than the densities of the massive subhalos in dark-matter- 
only simulations. SIDM offers a way to solve this problem be- 
cause it reduces the central density of halos. Thus in SIDM, the 
massive subhalos do host the luminous dSph but have shallower 
density profiles than predict ed in CDM simulations. T his has re- 
cently been demons trated bv | Vogelsberger et alj d2012l) . We do not 
directly compare to Vogelsberge r et al. d2012 ) because their work 
is focused on the subhalos of the Milky Way and the velocity- 
independent cross section that they simulate (a/m = 10cm 2 /g) 
is larger than the cross sections considered in our work. 

Regardless of whether Milky Way dSphs have cuspy or cored 
dark-matter halos, we may estimate the enclosed mass, and hence 
average density, around the half-light radius of the stellar dis- 
tribution. Mass estimates within 300pc and mass profile model- 
ings using stellar kinematics together with chemo-dynamically dis- 
tinct stellar subcomponets of Milky Way dwarf spheroidal galax- 
ies suggest central densities of app r oximately ~ 0.1Mp)/pc 3 
dStrigari et al.ll2008l: Iwolf et al.ll20ld IWalker & Pefiarrubidhoi ll; 
lAmorisco & Evansll2012l ; IWolf & Bullockll2012h . For the faintest 
dSph Segue 1, the density within the half-li ght radius (about 4 



pc) is measured to be about 2.5+ Mm /pc ISimon et al.1 d201 ll) ; 
iMartinez et alj d201 ll) . The errors on Segue 1 density are large but 
it is clear that if SIDM is to accommodate this result, it must al- 
low for large scatter in the core sizes and densities for small Vmxc 
halos. With a factor of 2-3 scatter in the densities quoted ear- 
lier for SIDMo.i halos, Segue 1 would appear to be compatible 
with SIDMo.i if its Vmax value is towards the lower end of the 
20 — 50 km/s range in V ma x. 

For the two dSph galaxies that appear to have cored 
density profiles (Fornax and Sculptor), the cores sizes must 
be of order- 0.2 - 1 kpc dWalker & Pefiarrubi j 1201 lh . For 
small halos with circular velocities in the 20 — 50 km/s range, 
which is close to the expected peak circular velocities of 
dwarf spheroidal halos before infall into the Milky Way host 
halo, an SIDM with a/m — 1 cm 2 /g predicts core sizes 
in the order of ~ 0.8 — 3.0 kpc, with central densities of 
about ~ 0.02 - 0.04M Q /pc 3 . Therefore, we find again that 
a/m = 1 cm 2 /g cannot reproduce the observed high central 
densities. On the other hand, our estimates suggest that an SIDM 
model with a/m — 0.1 cm 2 /g would produce central densities 
and core sizes consistent with the Milky Way dSph. 

In this last section we have used the analytic results that explain 
the scaling relations for the core sizes and central densities of halos 
in our SIDMi and SIDMo.i simulations, to extrapolate our results 
to scales ranging from galaxy clusters to dwarf spheroidal galaxies 
and to lower cross sections. We have found that a/m — 1 cm 2 /g 
would be unable to reproduce the observed high central densities. 
Remarkably, we find that the observations should be consistent with 
the predictions of a self-interacting dark matter with cross section 
in the ballpark of a/m — 0.1 cm 2 /g. These expectations are based 
on the scaling relations seen in SIDMi simulations and our ana- 
lytic model, which is consistent with the results from our direct 
a/m = 0.1 cm 2 /g simulations at the radii where we can trust our 
simulations. This deserves further study both in terms of simula- 
tions with SIDM cross section values smaller than 1 cm 2 /g and 
more detailed comparisons to observations. Our current look at the 
global data does not suggest a need for a velocity dependent cross 
section as has been previously suggested. In the companion paper 
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(Peter, Rocha, Bullock and Kaplinghat, 2012) we show that these 
SIDM models are also consistent with observations of halo shapes. 

8.3 Observed Substructure vs. SIDM 

In Figure [8] we showed that the number of subhalos for a/m = 
1 cm 2 /g is not significantly different from CDM predictions, 
especially in galaxy-scale halos. This is interesting because it 
means that SIDM fails to deliver on one of the original mo- 
tivations for conside r ing th is model of dark matter. Recall that 
ISpergel & Steinhardtl (2000) originally prom oted SIDM a s a so- 
lution to the m issing satellites problem dKlvpin et al.l 1 19991 ; 
Moor e et al .1 19991) . stating that many subhalos would be evaporated 
by interactions with the background halo. Given the new discover- 
ies of ultra-faint galaxies around the Milky Way and the hig h likeli- 
hood of many more discov eries from surveys like LSST dWillmanl 
l20ld : iBullocketalB oiO), a significant reduction in substructure 
counts may very well be a n egative characteristic of any non-CDM 
model dTollerud et alj2008l) . 

However, in Milky Way mass halos, SIDM with a/m = 
1 cm 2 /g will yield a significant probability for subhalo particle 
scattering only for systems that pass within ~ 10 kpc of the host 
halo center. Thus, for this cross section, we can form interesting- 
sized cores but largely leave the subhalo mass function unaffected 
in Milky Way-mass halos. For smaller cross sections, the differ- 
ences between SIDM and CDM subhalo mass functions will be 
even smaller. We note that we are not the first to find that SIDM 
can form cores b ut not solve the missing sate llites problem; it was 
first discussed in lD'Onghia & Burkertl d2003l) . 

This finding is also interesting in the context of other alter- 
natives to CDM. Warm dark matter (WDM) models, for which 
the outstanding difference from CDM is that dark-matter particles 
have high speeds at matter-radiation equality and a related free- 
streaming cutoff in the matter power spectrum, predict a suppres- 
sion in the halo (and subhalo) mass function at small scales. Oth- 
erwise, the abundance and stru cture of halos and subhalos is nearly 
indistinguishable fro m CDM dVillaescusa-Navarro & Dalai 1201 ll ; 
iMaccio' et alll2012h . WDM halos may be less concentrated than 
CDM halos on scales not much larger than the free-streaming scale, 
but are still cusped. They are only significantly cored right at the 
free-streaming scale, at which the halo and subhalo abundance is 
highly suppressed. Thus, each of the two leading modifications to 
CDM can solve only one of the two historical motivations for look- 
ing beyond the CDM paradigm. 

The lack of subhalo suppression for a/m < 1 cm 2 /g has im- 
plic ations for another of the S IDM halo constraints from a decade 
ago. lGnedin & O striker (2001) set a constraint excluding the range 
of 0.3 < a/m < 10 4 cm 2 /g based on the fundamental plane 
of elliptical galaxies. The argument rests on the observation that 
there are not significant differences in the fundamental plane of 
field ellipticals and cluster ellipticals (e.g., Kochanek et alj|200fl ; 
iBernardi et aT] 120031 ; lLa Barbera eTall l20ld) . Elliptical galaxies 
have a significant amount of dark matter within their half-light 
radii, with more massive ellipticals having larger mass-to-light ra- 
tios, either caused by varyi ng stellar mass-to-light ratios or vary- 
ing dark matter content jPa dmanabhan et al. 2004; Tolleru cTet al.l 
l201ll;IConrov & van Dokkumll2012l) . lGnedin & Ostrikedd200ll) "ar^ 
gue that elliptical galaxies falling into cluster-mass halos should 
have dark matter evaporated from their centers if a/m 7^ 0, which 
would cause the stars in the elliptical galaxy to adiabatically expand 
and hence move the galaxy off the fundamental plane. 

However, in our simulations, we find that few subhalos are 



fully evaporated, and that the subhalo Vm ax function is not greatly 
different for a/m — 1 cm 2 /g from CDM. In addition, our ana- 
lytic arguments show that the trend with (host) halo mass for the 
eva poration of subhalos at fix ed r/r s is mild. This suggests that 
the iGnedin & OstrikeJ d200ll) constraints are overly conservative 
even at the a/m ~ 1 cm 2 /g level. The main caveats are that the 
suppression of the subhalo V mB function is higher in more mas- 
sive clusters and that the suppression is highest at the center of 
the cluster halo. It would also be interesting to see if there are any 
differences in the fundamental plane as a function of projected dis- 
tance in the cluster, both observationally and in simulations. For 
all of these reasons, it would be worthwhile to perform simula- 
tions of elliptical galaxies in clusters with SIDM and explore the 
fundamental-plane constraints in more depth. 

To summarize, although we have not fully resolved the 
cores of a/m = 0.1cm 2 /g SIDM halos, the intuition gleaned 
from our analytic model (tested agains the SIDMi results) and 
our moderately-resolved simulation results suggest that a/m — 
0.1 cm 2 /g is an excellent fit to the data across the range of halo 
masses from dwarf satellites of the Milky Way to clusters of galax- 
ies. Values of cross section over dark matter particle mass in this 
range are fully consistent with the published Bullet cluster con- 
straints (cf. fT), measurements of dark matter density on small- 
scales and subhalo survival requirements. In a companion paper 
(Peter, Rocha, Bullock and Kaplinghat 2012), we show that this 
model is also consistent with halo shape estimates. It is therefore 
important to simulate galaxy and cluster halos with cross sections 
in the 0.1 cm 2 /g range. 



9 SUMMARY AND CONCLUSIONS 

We have presented a new algorithm to include elastic self-scattering 
of dark matter particles in N-body codes and used it to study the 
structure of self-interacting dark matter (SIDM) halos simulated 
in a full cosmological context. Our suite of simulations (summa- 
rized in Table 1) rely on identical initial conditions to explore SIDM 
models with velocity-independent cross sections a/m = 1 cm 2 /g 
and a/m — 0.1 cm 2 /g as well as a comparison set of standard 
CDM simulations (with a/m = 0). 

Our primary conclusion is that while SIDM looks identical 
to CDM on large scales, SIDM halos have constant density cores, 
with core radii that scale in proportion to the standard CDM scale 
radius (r corc ~ er s ). The relative size of the core increases with 
increasing cross section (e ~ 0.7 for a/m = 1 and e ~ 0.2 for 
a/m — 0.1 cm 2 /g). Correspondingly, at fixed halo mass, core 
densities decrease with increasing SIDM cross section. For both 
core radii and core densities, there is significant scatter about the 
scaling with V max of the halo. The scaling relationship is strong 
enough that measurements of dark matter densities in the cores 
of dark matter dominated galaxies and large galaxy clusters likely 
provide the most robust constraints on the dark matter cross sec- 
tion at this time. In a companion paper (Peter, Rocha, Bullock and 
Kaplinghat, 2012) we demonstrate, contrary to previous claims, 
that SIDM constraints from halo shape measurements may be less 
restrictive than (or at least similar to those from) measurements of 
absolute core densities alone. 

Based on our simulation results we conclude that the dark 
matter self-scattering cross section must be smaller than 1 cm 2 /g 
in order to avoid under-predicting the observed core densities in 
galaxy clusters, low surface brightness spirals (LSBs), and dwarf 
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spheroidal galaxies. However, an SIDM model with a velocity- 
independent cross section of about a/m = 0.1 cm 2 /g appears 
capable of reproducing reported core sizes and central densities of 
dwarfs, LSBs, and galaxy clusters. Higher resolution simulations 
with better statistics will be needed to confirm this expectation. 

An accounting of our results are as follows: 

• Outside of the central regions of dark matter halos 
(r > 0.5i?vir) the large scale properties of SIDM cosmolog- 
ical simulations are effectively identical to CDM simulations. This 
implies that all of the large-scale confirmations of the CDM theory 
apply to SIDM as well. 

• The subhalo V max function in SIDM with a/m — 1 cm 2 /g 
differs by less than ~ 30% compared to CDM across the mass 
range 5 x 1O 11 M — 2 x 1O 14 M studied directly with our simu- 
lations . Differences in the V max function with respect to CDM are 
only apparent deep within the centers of large dark- matter halos. 
Thus, although is possible, it will be difficult to constrain SIDM 
models based on the effects subhalo evaporation. 

• SIDM produces halos with constant density cores, with 
correspondingly lower central densities than CDM halos of the 
same mass. For a/m = 1 cm 2 /g, our simulated halo density 
structure is reasonably well characterized by a Burkert (1995) 
profile fit with a core size ?~b — 0.7r s , where r s is the NFW scale 
radius of the same halo in the absence of self-interactions. Core 
densities tend to increase with decreasing halo mass (pb oc M~® 2 ) 
but demonstrate about a factor of ~ 3 scatter at fixed mass (likely 
owing to the intrinsic scatter in dark matter halo concentrations). 

• SIDM halo core sizes, central densities, and associated 
scaling relations can be understood in the context of a simple 
analytic model. The model treats the SIDM halo as consisting of 
a core region, where self-interactions have redistributed kinetic 
energy to create an approximately isothermal cored density profile; 
and an outer region, where self-interactions are not effective. 
The transition between these regions is set by the strength of the 
self-interactions and this model allows us to make quantitative 
predictions for smaller cross sections where the cores are not 
resolved by our simulations. Based on this model and a few of 
our best resolved simulated halos we find core sizes ~ 0.1r s for 
a/m — 0.1 cm 2 /g. 

• Halo core densities over the mass range from 
1Q i5 _ 10 io Mq in SIDM with a j m = i cm 2 /g are too 

low (~ 0.005 — 0.04 Mq/pc 3 ) to match observed central den- 
sities in galaxy clusters (~ 0.03M Q /pc) and dwarf spheroidals 
(~O.lM /pc 3 ). 

• Halo core central densities in SIDM with a/m = 0.1 cm 2 /g 
are in line with those observed from galaxy clusters to tiny dwarfs 
(0.02 — 0.5M Q /pc 3 ) without the need for any velocity depen- 
dence. The densities are more consistent with observations than 
those predicted in dissipationless CDM simulations, which are 
generically too high. SIDM models with this cross section over 
dark matter particle mass value are consistent with Bullet cluster 
observations, subhalo survival requirements and, as we show in a 
companion paper (Peter, Rocha, Bullock and Kaplinghat, 2012), 
measurements of dark matter halo shapes. 

Future work is necessary to expand both the dynamic range of 



our simulations in halo mass and resolution as well as the dynamic 
range in cross sections. These simulations are necessary in order 
to make detailed comparisons with observations given the exciting 
possibility that dark matter self-interaction with a/m in the ball- 
park of 0.1 cm 2 /g could be an excellent fit to the central densities 
of halos over 4-5 orders of magnitude in mass. 
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APPENDIX A: DERIVATION OF THE HARD-SPHERE INTERACTION RATE IN N-BODY SIMULATIONS 

The challenge is to represent a microphysical scattering process in a macroscopic context in which neither a fluid nor collisionless treatment 
is appropriate. In order to develop a Lagrangian technique in which to represent the scattering process, we start with the Boltzmann equation. 
Particles with mass m, a hard-sphere scattering cross section da/dQ (as a function of center-of-mass scattering angle), and a distribution 
function /(x, v, t) evolve as 

D/(x,v,t) =lUff] (A1) 



Dt 



J d 3 Vl J dfi^lv-vrl [/(x,v',i)/(x,vi,t)-/(x,v,t)/(x,v 1 ,i)] . (A2) 



Here, D/Dt is a Lagrangian time derivative and T[f, a] is the collision operator. If the particles were collisionless, the Lagrangian time- 
derivative of the distribution function would be zero; the phase-space density of particles would be conserved. The left-hand expression in 
the brackets in Equation dA2t represents scattering of particles into a small patch of phase space centered on (x, v), and the right-hand 
expression (after the — sign) represents scattering out of that patch of phase space. If v and vi represent the initial velocities of the primary 
and target particles, then v' and are their post-scatter velocities, which are related to the initial velocities by the center-of-mass scattering 
angle SI. 

The key step in being able to represent the scattering process in a simulation is the ansatz that the evolution of the coarse-grained 
distribution function / (the distribution function averaged over several times the interparticle spacing) is a good representation of the evolution 
of the fine-grained distribution function /. In other words, the ansatz is that the solution to 



5* = / ^ / d "S |v_Vl1 [A«.v',t)/(x,vi,t)-/(x,v,t)/(x,vi,t)] (A3) 



is the same as the solution for / in Equation dA2t averaged over a patch of phase space. If this is the case, our next step is to discretize 
Equation dA3b such that we can solve the Boltzmann equation by Monte Carlo N-body methods. 

To discretize Equation dA3l >. we consider a particle-based Lagrangian method in which each particle in the N-body simulation represents 
a patch of phase space. In the absence of collisions, the simulation particles trace out geodesies in the gravitational field of the particles. When 
we discretize the phase space, we do it as follows: 

/(x,v,i) = ^(M i /m)W(|x-x i |;/i i )5 3 (v-v i ). (A4) 

i 

Here, i labels a discrete macro particle representing a patch of phase space that has mass Ah ; thus each macro particle represents a patch 
of phase space inhabited by Mi/m of the true particles. We assume a delta-function form for the velocity distribution because each macro 
particle travels at only one speed. We treat each macro particle as being smoothed out in configuration space with a smoothing kernel W 
with smoothing length hi . The reason for treating each macro particle as inhabiting a finite region of configuration space is that we want the 
local estimate of the density 



n{x) = J 



d 3 v/ (A5) 



to be smooth. Preliminary tests show that smoothness is necessary to properly estimate the collision term of the Boltzmann equation. Note 
that in the main text we use Mi = m p and hi = h s i for all i. This is because all of the N-body particles have the same mass in our simulations 
and we have fixed h s i to be constant for all particles in the simulations we present. 

In the particle-based discretization of the Boltzmann equation, the fact that each particle represents a patch of phase space means that 
we must discretize the collision operator; we must integrate the collision operator over the patch of phase space inhabited by a single particle. 
Thus, if a specific particle represents a patch of phase space of size 5x p 5v p , we must calculate 



j d 3 x^ ^%=f d 3 x^ d B vyd 8 viy'dn^|v-vi|[/(x,v , ,t)/(x,vi,t)-/(x,v,t)/(x,vi,t)] 



(A6) 



Thus, our approach to estimating the collision term and the Boltzmann equation are as follows. To find the collision rate for the region of 
phase space associated with a particle j, we divide Equation dA6t through by Mj/m (so that we are calculating the scattering probability 
for a single macro particle j), and we consider only the "scattering out" part of the collision operator. We consider the pairwise rate Tij for 
particle j to scatter off any of the other i particles. We do a Monte Carlo simulation of the scatters; if a pair of particles is allowed to scatter 
in a given small timestep, we calculate the macro particles' post-scatter velocity using the center-of-mass scattering angle £1. This latter step 
is our approximation to the "scatter out" term of the Boltzmann collision operator. 
The pairwise collision operator is 

_ r( P \q) + r( q \p) 

1 pq — o 1 



where the conditional probability of scattering a specific particle p off a target particle q is T(q\p), which is determined by the collision term 
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of the Boltzmann equation. This collision term is derived from Equation dA6l l. such that 

r(p)= f d 3 x f d\ f d 3 Vl f dn^\v~v 1 \(M p /m)- 1 f(x,v,t)f(x,v 1 ,t) (A8) 

J Sxp J Svp J J 

= f d 3 x / d 3 v / d 3 vi /dfi^|v-vi|(Mj,/m)" 1 {V(A/ j /m)W r (|x-x j |;/i j )<5 3 (v-v J ) 
iix p iiv p J J all - 

^2{M q /m)W(\x - x,|; /i,)5 3 ( Vl - v q )} (A9) 

= ^ d 3 x y d 3 Vl y'dn^|v p -v 1 |^(M,/m) 1 ^(!x-x p |;/i p )W(|x-x 9 |;/i 9 )5 3 (v 1 -v q ) (A10) 

= E f dQ %^r^«-^ I rf 3 xVK(|x-x p |;/ tp )^(|x-x 9 |;^) (All) 

= J2( a / m ) M i\ v i - v p\9pi ( A12 > 

9 

= ^r( g |p). (Ai3) 

9 

We note the appearance of the term a/m, which is the scattering cross section per unit mass. The kernel g is defined as 

pmax (hp ,hq ) 

g Pq = d 3 x'W(\x'\,h p )W(\Sx pq + x'\,h q ). (A14) 

Jo 

Using these sets of equations, we calculate T pq for each pair of particles whose configuration-space patches overlap at each time step St, 
making sure to keep the time steps small enough that T pq 5t <C 1 for each time step. 

APPENDIX B: TEST FOR THE SCATTERING KINEMATICS 

We use the same setup described in ^3]to test our implementation against the expected kinematics. For this we looked at the distributions of the 
post-scatter velocity magnitudes and directions for both the Sphere and Background particles. For the distributions of the velocity directions 
we looked at the inclination and azimuthal angles of the post-scatter velocity vectors. The angles are defined such as the line of interaction is 
along the 9 — direction (i.e the z-axis) and <j> is the azimuthal angle about which the experiment is symmetric. The distributions resulting 
from our test simulation are compared to those obtained from the transformation of a uniform isotropic distribution in the center of mass 
frame to the simulation/lab frame; the results are shown in Figures lBllB3l Figure lBll shows that the distributions of the velocity magnitudes 
rise linearly from v = to v — v s , followed by a sharp cutoff at v = v s , where v s is the relative speed between the sphere and the 
background. From conservation of energy it is only possible to have particles with v > v 3 if they have interacted multiple times. Multiple 
interactions are not considered in our calculations of the theoretical distributions but they are allowed in our simulation, hence, in Figure 
IB II one can observe a tail for velocities > v s on the distributions of both types of particle velocities, but not on the theoretical distribution. 
Looking at Figure lB2l one can see that most of the particles are scattered towards the 9 — 45° directions, i.e. forming a 45° angle with respect 
to v s . It is visible that the distributions resulting from the simulation in the left panel of Fi sure lB2l are higher than expected for 9 <~ 20°. 
This is again from the fact that multiple scatters are possible in the simulation and this is not included in the calculations of the theoretical 
histograms. We demonstrate that this is the case by showing in the right panel of Figure |B2] the distributions obtained from the simulation 
when we exclude any particles with v > v s , excluding that way any particles that we know have interacted multiple times, as one can see 
doing this brings the distributions from the simulation to a better match with the theory. Figure |B3] shows the distributions of the velocities 
as a function of (f>, these are flat as expected due to the symmetry of the experiment. 
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Figure Bl. Distribution of the post-scatter velocity magnitudes. From conservation of energy it is only possible to have particles with velocities > v s if they 
have interacted multiple times, this is not included in our calculation of the theoretical distribution but it is allowed in our simulation, hence one can observe a 
tail for velocities > v 3 on the distributions of both types of particle velocities, but not on the theoretical distribution. 




Figure B2. Distributions of the post-scatter velocities along the ^-directions. It is evident that most of the particles are scattered towards the 9 = 45° directions, 
i.e. forming a 45° angle with v s . Note that the distributions resulting from the simulation in the left panel are higher that expected for 8 <~ 20°. This is 
because multiple scatters are possible in the simulation and they are not considered in the calculations of the theoretical histograms. We demonstrate this by 
showing in the right panel the distributions from the simulation when we exclude any particles with v > v s , excluding that way any particles that we know 
have interacted multiple times and bringing the distributions from the simulation to a better agreement with the theory. 
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Figure B3. Distributions of the velocities along the (^-directions. The flat distributions show that the results are symmetric about the direction of motion, i.e. 
the z-axis. 



